|
| 1 | +--- |
| 2 | +title: Hail GVCF Combiner Standard Operating Procedure |
| 3 | +author: Chris Vittal |
| 4 | +date: 2026-09-02 |
| 5 | +updated: 2026-09-02 |
| 6 | +--- |
| 7 | + |
| 8 | +# Hail GVCF/VDS Combiner Standard Operating Procedure |
| 9 | + |
| 10 | +## Contents |
| 11 | + |
| 12 | +1. [Overview](#overview) |
| 13 | +1. [Collating Samples](#collating-samples) |
| 14 | +1. [Ensuring Access](#ensuring-access) |
| 15 | +1. [Combining](#combining) |
| 16 | +1. [Troubleshooting](#troubleshooting) |
| 17 | + |
| 18 | +## Overview |
| 19 | + |
| 20 | +This document is intended to describe a procedure for the creation of a hail |
| 21 | +Variant Data Set (VDS) from a collection of input GVCFs. It details best |
| 22 | +practices and gives examples for using Hail Batch on GCP to perform this |
| 23 | +work. |
| 24 | + |
| 25 | +## Collating Samples |
| 26 | + |
| 27 | +In order to create a joint dataset, there must first be single sample data. This |
| 28 | +will usually be provided by a lab project manager and contain for every sample, |
| 29 | +at minimum, the desired final sample ID (more on this in a bit), and the object |
| 30 | +storage URI for that sample. |
| 31 | + |
| 32 | +It is an exercise for the operator to take that data and turn it into two python |
| 33 | +lists of the same length, one a list of sample IDs and the other the list of |
| 34 | +object URIs. |
| 35 | + |
| 36 | +### The `new_combiner` API (part 1) |
| 37 | + |
| 38 | +Before going further, it is useful to begin understanding the inputs to |
| 39 | +`new_combiner`, the public constructor for the `VariantDatasetCombiner` |
| 40 | +class (the combiner driver). |
| 41 | + |
| 42 | +```python3 |
| 43 | +def new_combiner( |
| 44 | + *, |
| 45 | + gvcf_paths: Optional[List[str]] = None, |
| 46 | + gvcf_sample_names: Optional[List[str]] = None, |
| 47 | + vds_paths: Optional[List[str]] = None, |
| 48 | + vds_sample_counts: Optional[List[int]] = None, |
| 49 | + ..., |
| 50 | +) -> VariantDatasetCombiner |
| 51 | +``` |
| 52 | + |
| 53 | +Above are the four parameters that describe the inputs to the combiner. Of |
| 54 | +these, only one of `gvcf_paths` or `vds_paths` is actually required. Operators |
| 55 | +of the combiner however SHOULD when providing `gvcf_paths`, provide |
| 56 | +`gvcf_sample_names` and when providing `vds_paths`, provide `vds_sample_counts`. |
| 57 | +This will save work on the backend as the combiner pipelines compute the sample |
| 58 | +counts and sample names if necessary. |
| 59 | + |
| 60 | +`gvcf_sample_names` and `gvcf_external_header` must be set together, or both |
| 61 | +left unset. Setting only one raises a `ValueError`. See |
| 62 | +[Controlling Imported Fields](#controlling-imported-fields) for `gvcf_external_header`. |
| 63 | + |
| 64 | +This is particularly important for GVCFs as often two samples across projects or |
| 65 | +workspaces can have the same ID (and therefore there will be duplicate sample |
| 66 | +IDs in the final VDS), but should be distinguished usually with some kind of |
| 67 | +project prefix. |
| 68 | + |
| 69 | +## Ensuring Access |
| 70 | + |
| 71 | +It is _annoyingly_ common to waste time due to cloud permissions errors. Because |
| 72 | +of this, operators should check every input—both vcf and tabix—to ensure they |
| 73 | +have access to all of them. |
| 74 | + |
| 75 | +The following python function should suffice to check a list of URIs. When using |
| 76 | +it, remember to include the `.tbi` tabix index file for each input GVCF. |
| 77 | +A comprehension like `[file for gvcf in gvcfs for file in (gvcf, gvcf + '.tbi')]` |
| 78 | +will add the corresponding index for each file. |
| 79 | + |
| 80 | +```python3 |
| 81 | +import sys |
| 82 | + |
| 83 | +from hailtop.aiotools.fs import AsyncFS |
| 84 | +from hailtop.utils import bounded_gather |
| 85 | + |
| 86 | +async def check_permissions(fs: AsyncFS, files: list[str]): |
| 87 | + result = await bounded_gather(*(lambda f=file: fs.statfile(f) for file in files), |
| 88 | + parallelism=50, |
| 89 | + return_exceptions=True) |
| 90 | + missing_permissions = [(file, exc) for file, (_, exc) in zip(files, result) if exc is not None] |
| 91 | + if not missing_permissions: # no exceptions, every file can be accessed |
| 92 | + return |
| 93 | + |
| 94 | + for file, exc in missing_permissions: |
| 95 | + print(f"error reading {file}, exception =", exc, file=sys.stderr) |
| 96 | + msg = f"Missing permissions (or files don't exist) for {len(missing_permissions)} of {len(files)} files" |
| 97 | + raise PermissionError(msg) from missing_permissions[0][1] |
| 98 | +``` |
| 99 | + |
| 100 | +## Combining |
| 101 | +Once all GVCFs have been collated and permissions checked it is finally time to |
| 102 | +run the combine job. |
| 103 | + |
| 104 | +Actually running the combiner is as simple as: |
| 105 | +```python3 |
| 106 | +combiner = hl.vds.new_combiner(...) |
| 107 | +combiner.run() |
| 108 | +``` |
| 109 | + |
| 110 | +To actually create a combiner object, we walk through each remaining portion of |
| 111 | +the `new_combiner` constructor giving good defaults and recommendations. |
| 112 | + |
| 113 | +### Required and Recommended Parameters |
| 114 | + |
| 115 | +```python3 |
| 116 | +def new_combiner( |
| 117 | + *, |
| 118 | + output_path: str, |
| 119 | + temp_path: str, |
| 120 | + save_path: Optional[str] = None, |
| 121 | + ..., |
| 122 | +) -> VariantDatasetCombiner: |
| 123 | +``` |
| 124 | + |
| 125 | +The only two specifically required parameters to `new_combiner` are `output_path` |
| 126 | +and `temp_path`. The operator MUST set `output_path` to a durable object storage |
| 127 | +URI (no lifecycle policy on a bucket). The operator MUST set `temp_path` to |
| 128 | +a transient object storage URI (a bucket with a lifecycle policy), however it is |
| 129 | +RECOMMENDED that the lifecycle policy be long enough to keep all created objects |
| 130 | +live for the duration of the whole combine job—usually up to two weeks. This |
| 131 | +ensures that intermediates are always available while also ensuring timely |
| 132 | +cleanup to keep storage costs down. |
| 133 | + |
| 134 | +The operator SHOULD set `save_path` to an object storage URI, this can be either |
| 135 | +temporary or durable. It is RECOMMENDED that it be temporary. This parameter |
| 136 | +controls wher the combiner object serializes itself as JSON. This can and will |
| 137 | +be loaded to pick up where combining left off in the event of an error like an |
| 138 | +OOM or cloud error. The combiner will produce a path based on `temp_path` and |
| 139 | +the hash of its parameters. However, this is liable to change without notice so it |
| 140 | +is best to set a path for easy introspection on the plan and ease of use. |
| 141 | + |
| 142 | +### Import VCF parameters |
| 143 | + |
| 144 | +```python3 |
| 145 | +def new_combiner( |
| 146 | + *, |
| 147 | + reference_genome: Union[str, ReferenceGenome] = 'default', |
| 148 | + contig_recoding: Optional[Dict[str, str]] = None, |
| 149 | + call_fields: Collection[str] = ['PGT'], |
| 150 | + ..., |
| 151 | +) -> VariantDatasetCombiner: |
| 152 | +``` |
| 153 | + |
| 154 | +The operator SHOULD set `reference_genome` to the appropriate reference genome |
| 155 | +to avoid issues where the default does not match the reference of the samples. |
| 156 | + |
| 157 | +To understand `contig_recoding`, refer to the [`import_vcf`] documentation. The |
| 158 | +operator MUST set this parameter if contigs that do not match hail's |
| 159 | +representation of the reference genome have been used for the GVCFs. |
| 160 | + |
| 161 | +Use `call_fields` if there are fields other than `GT` and `PGT` that should be |
| 162 | +treated as a genomic call, again refer to [`import_vcf`] for more details. |
| 163 | + |
| 164 | +[`import_vcf`]: https://hail.is/docs/0.2/methods/impex.html#hail.methods.import_vcf |
| 165 | + |
| 166 | +### Interval Parameters |
| 167 | + |
| 168 | +```python3 |
| 169 | +def new_combiner( |
| 170 | + *, |
| 171 | + intervals: Optional[List[Interval]] = None, |
| 172 | + import_interval_size: Optional[int] = None, |
| 173 | + use_genome_default_intervals: bool = False, |
| 174 | + use_exome_default_intervals: bool = False, |
| 175 | + ..., |
| 176 | +``` |
| 177 | + |
| 178 | +When combining GVCFs (`gvcf_paths` is a nonempty list), exactly one of the four |
| 179 | +interval parameters MUST be set. Though not prohibited by the code it is still |
| 180 | +incorrect and can lead to surprising behavior (a warning will be raised |
| 181 | +otherwise, and the parameter first in the above set of arguments that is set |
| 182 | +will be used). It is usually sufficient to use one of the |
| 183 | +`use_*_default_intervals` parameters. These give reasonably even coverage over |
| 184 | +the genome for each sample type (exome/genome). When combining VDS files only, |
| 185 | +the interval parameters are ignored and `intervals` is set to an empty list. |
| 186 | + |
| 187 | +The point of the interval parameters is to control parallelism. The default |
| 188 | +exome intervals are rather coarse since exomes are small (around 70). The |
| 189 | +default genome intervals include a few thousand items. |
| 190 | + |
| 191 | +When setting `intervals`, each interval provided MUST be non-overlapping, on |
| 192 | +a single contig, and include both endpoints. For a production run, `intervals` |
| 193 | +SHOULD cover almost all of the canonical contigs for the samples' reference |
| 194 | +genome (excluding some amount of telomeres/centromeres is often ok). The operator |
| 195 | +MUST NOT use "calling intervals" to provide `intervals`, exome calling intervals |
| 196 | +are much finer grained and will lead to far too much parallelism. Use |
| 197 | +`intervals` when large regions are to be excluded, or for testing as the |
| 198 | +explicit `intervals` parameter is the only one that can filter out large |
| 199 | +regions. For example, when testing, only intervals from chr20 can be used to |
| 200 | +avoid reading all of every sample. |
| 201 | + |
| 202 | +Setting `import_interval_size`, causes the library to compute even intervals of |
| 203 | +roughly that size to then use for GVCF import. Use this if a default is too |
| 204 | +coarse or fine for the samples involved. (It is difficult to define 'too much' |
| 205 | +or 'too little' parallelism, if tasks are taking too long, more |
| 206 | +intervals/smaller interval size can be used; if there is too much overhead due |
| 207 | +to excess parallelism a larger interval size can be used.) |
| 208 | + |
| 209 | +### Controlling Imported Fields |
| 210 | + |
| 211 | +There are four parameters that help control what fields are kept/discarded on |
| 212 | +GVCF import. |
| 213 | + |
| 214 | +```python3 |
| 215 | +def new_combiner( |
| 216 | + *, |
| 217 | + gvcf_external_header: Optional[str] = None, |
| 218 | + gvcf_info_to_keep: Optional[Collection[str]] = None, |
| 219 | + gvcf_reference_entry_fields_to_keep: Optional[Collection[str]] = None, |
| 220 | + gvcf_save_filters: bool = False, |
| 221 | + ..., |
| 222 | +) -> VariantDatasetCombiner: |
| 223 | +``` |
| 224 | + |
| 225 | +It is RECOMMENDED that the operator set `gvcf_external_header`. This serves two |
| 226 | +purposes. The first is that it instantly homogenizes the input type of each gvcf |
| 227 | +thus avoiding platform/software differences in gvcf creation. Second, without |
| 228 | +this parameter, every file needs to be checked to ensure homogeneity, so setting |
| 229 | +it saves some time. This will not cause issues, due to the way that VCF |
| 230 | +parsing in hail works, fields that are not in the header, but are present in the |
| 231 | +VCF are silently skipped during parsing. Before setting this parameter, confirm |
| 232 | +that the header declares every `INFO` and `FORMAT` field the analysis needs. |
| 233 | + |
| 234 | +The operator SHOULD set `gvcf_save_filters` to `True` in order to save any |
| 235 | +filters present in the GVCF; it defaults to `False` for backwards compatibility. |
| 236 | + |
| 237 | +It is not usually necessary to set `gvcf_info_to_keep` or |
| 238 | +`gvcf_reference_entry_fields_to_keep`, the combiner does a good job of |
| 239 | +determining which fields are present in reference/variant data and keeping them. |
| 240 | +There is one notable field that operators may wish to keep in reference data: |
| 241 | +`PL` is dropped by default for reference data and must be explicitly added to |
| 242 | +`gvcf_reference_entry_fields_to_keep` if it is desired to be retained. |
| 243 | +Unfortunately setting `gvcf_reference_entry_fields_to_keep=['PL']` will drop |
| 244 | +everything _but_ `PL`, so the other fields to be kept must be added manually as |
| 245 | +well. |
| 246 | + |
| 247 | +### Parallelism Control |
| 248 | +```python3 |
| 249 | +def new_combiner( |
| 250 | + *, |
| 251 | + branch_factor: int = 100, |
| 252 | + target_records: int = 24_000, |
| 253 | + gvcf_batch_size: int = 50, |
| 254 | + ..., |
| 255 | +) -> VariantDatasetCombiner: |
| 256 | +``` |
| 257 | + |
| 258 | +These parameters work to control how many inputs get combined at once |
| 259 | +(`branch_factor`, default 100), how large final output partitions should be |
| 260 | +(`target_records`, default 24,000), and how many gvcf combine jobs to attempt to |
| 261 | +run in parallel (`gvcf_batch_size`, default 50). |
| 262 | + |
| 263 | +These are the only parameters that do not get taken into account when saving the |
| 264 | +combiner plan to `save_path`. This is so that they can be changed mid-run, as |
| 265 | +tuning these values is often essential to making sure the combiner finishes—more |
| 266 | +on that later. |
| 267 | + |
| 268 | +Operators SHOULD leave `branch_factor` and `gvcf_batch_size` as defaults and |
| 269 | +adjust them as part of troubleshooting. It is RECOMMENDED that operators work |
| 270 | +with downstream analysts to understand parallelism needs and set target records |
| 271 | +accordingly, values up to a few million are not inappropriate for lower numbers |
| 272 | +of samples (less than 10k) in a VDS. |
| 273 | + |
| 274 | +## Troubleshooting |
| 275 | +The combiner is often a long running pipeline. It is resilient but may need to |
| 276 | +be tweaked or restarted. |
| 277 | + |
| 278 | +### Resources |
| 279 | +Combiner pipelines are very large, often using lots of memory for compilation. |
| 280 | +The jobs can fail with out of memory errors. Many strange errors of unknown |
| 281 | +cause are resource exhaustion errors in disguise. Increasing `driver_memory` or |
| 282 | +`worker_memory` in `hl.init` may be necessary for the jobs to proceed. |
| 283 | + |
| 284 | +Alternatively, increasing `target_records` can decrease memory used during |
| 285 | +compilation and execution by decreasing the number of intervals that each |
| 286 | +intermediate VDS is read with and therefore also the size of the array of |
| 287 | +partition results. |
| 288 | + |
| 289 | +### Retry |
| 290 | +The Hail Batch service has great preemption and retry logic, however it still |
| 291 | +can't catch everything. If there is a service error or a spurious/rare cloud |
| 292 | +error it is usually enough to restart the combine job. To do so, run the same |
| 293 | +combine script, it will load the combiner plan and continue. |
| 294 | + |
| 295 | +### Compilation errors |
| 296 | +The combiner generally only sees one form of compilation error, |
| 297 | +`ClassTooLargeException`. This usually occurs due to too many inputs and so is |
| 298 | +most often resolved by shrinking the number of inputs at the cost of |
| 299 | +parallelism. If this error is encountered in GCVF combining, try decreasing |
| 300 | +`gvcf_batch_size` by 25% until the pipeline runs as expected. If there are |
| 301 | +issues combining VDS files, try increasing `target_records` to start, especially |
| 302 | +if there are hundreds of thousands of partitions in the VDS combine jobs. Try |
| 303 | +doubling it to see if that makes a difference. |
| 304 | + |
| 305 | +As a last resort to keep the combiner working, `branch_factor` can be decreased. |
| 306 | +This parameter determines how many times the whole dataset needs to be read from |
| 307 | +start to finish so should be kept as large as feasible. Like batch size, try |
| 308 | +decreasing `branch_factor` by 25% until the pipeline runs again. |
0 commit comments