Repository: whatshap/whatshap

Issue:
whatshap phase: Make --indels the default
I suggest to enable indel phasing by default. The option doesn’t have the best name anyway (it also works with "complex" variants/MNPs). By instead adding an `--only-snvs` option, that naming issue would also be solved.

This may require bumping the version to 2.0.

Discussion:

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 276bc1e7..7cb126cf 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -5,6 +5,16 @@ Changes
 development version
 -------------------
 
+* :issue:`346`: Phasing of indels (and other non-SNVs) is now enabled by
+  default. This previously required specifying the ``--indels`` option,
+  which not all users knew about and were thus unnecessarily getting
+  suboptimal phasing results. The option is now ignored and leads to a
+  warning. An ``--only-snvs`` option was added that restores the old behavior.
+  This change applies to the following subcommands: ``phase``, ``haplotype``,
+  ``polyphase``, ``polyphasegenetic``.
+
+  Since this is a backwards incompatible change (when not using ``--indels``
+  already), the major version has been increased.
 * :issue:`425`: Haplotagging CRAM files should now work in more cases with
   ``haplotag``.
 * :issue:`457`: ``haplotag`` now also tags alignments marked as duplicate.
diff --git a/whatshap/cli/__init__.py b/whatshap/cli/__init__.py
index 3f26928a..2ae7498a 100644
--- a/whatshap/cli/__init__.py
+++ b/whatshap/cli/__init__.py
@@ -46,7 +46,7 @@ class PhasedInputReader:
         reference,
         numeric_sample_ids,
         ignore_read_groups,
-        indels,
+        only_snvs,
         **kwargs,  # passed to ReadSetReader constructor
     ):
         self._bam_paths, self._vcf_paths = self._split_input_file_list(bam_or_vcf_paths)
@@ -55,7 +55,7 @@ class PhasedInputReader:
         self._numeric_sample_ids = numeric_sample_ids
         self._fasta = self._open_reference(reference) if reference else None
 
-        vcf_readers = [VcfReader(f, indels=indels, phases=True) for f in self._vcf_paths]
+        vcf_readers = [VcfReader(f, only_snvs=only_snvs, phases=True) for f in self._vcf_paths]
 
         self._vcf_readers = vcf_readers
         self._ignore_read_groups = ignore_read_groups
diff --git a/whatshap/cli/compare.py b/whatshap/cli/compare.py
index 386c03d7..5b6c789e 100644
--- a/whatshap/cli/compare.py
+++ b/whatshap/cli/compare.py
@@ -763,7 +763,7 @@ def run_compare(
     longest_block_tsv=None,
 ):
     vcf_readers = [
-        VcfReader(f, indels=not only_snvs, phases=True, ploidy=ploidy, mav=True) for f in vcf
+        VcfReader(f, only_snvs=only_snvs, phases=True, ploidy=ploidy, mav=True) for f in vcf
     ]
     if names:
         dataset_names = names.split(",")
diff --git a/whatshap/cli/genotype.py b/whatshap/cli/genotype.py
index 33ce0f7a..f08f89ef 100644
--- a/whatshap/cli/genotype.py
+++ b/whatshap/cli/genotype.py
@@ -7,7 +7,8 @@ forward backward algorithm.
 import logging
 import sys
 import platform
-from typing import Sequence
+from argparse import SUPPRESS
+from typing import Sequence, Optional
 
 from contextlib import ExitStack
 from whatshap import __version__
@@ -74,7 +75,7 @@ def run_genotype(
     samples=None,
     chromosomes=None,
     ignore_read_groups=False,
-    indels=True,
+    only_snvs=False,
     mapping_quality=20,
     max_coverage=15,
     nopriors=False,
@@ -115,7 +116,7 @@ def run_genotype(
                 reference,
                 numeric_sample_ids,
                 ignore_read_groups,
-                indels=indels,
+                only_snvs=only_snvs,
                 mapq_threshold=mapping_quality,
                 overhang=overhang,
                 affine=affine_gap,
@@ -131,7 +132,7 @@ def run_genotype(
             GenotypeVcfWriter(command_line=command_line, in_path=variant_file, out_file=output)
         )
         # vcf writer for only the prior likelihoods (if output is desired)
-        prior_vcf_writer = None
+        prior_vcf_writer: Optional[GenotypeVcfWriter] = None
         if prioroutput is not None:
             prior_vcf_writer = stack.enter_context(
                 GenotypeVcfWriter(
@@ -145,7 +146,7 @@ def run_genotype(
         # remove all likelihoods that may already be present
         vcf_reader = stack.enter_context(
             VcfReader(
-                variant_file, indels=indels, genotype_likelihoods=False, ignore_genotypes=True
+                variant_file, only_snvs=only_snvs, genotype_likelihoods=False, ignore_genotypes=True
             )
         )
 
@@ -261,7 +262,7 @@ def run_genotype(
 
             # if desired, output the priors in separate vcf
             if prioroutput is not None:
-                prior_vcf_writer.write_genotypes(chromosome, variant_table, indels)
+                prior_vcf_writer.write_genotypes(chromosome, variant_table, only_snvs)
 
             # Iterate over all families to process, i.e. a separate DP table is created
             # for each family.
@@ -367,7 +368,7 @@ def run_genotype(
 
             with timers("write_vcf"):
                 logger.info("======== Writing VCF")
-                vcf_writer.write_genotypes(chromosome, variant_table, indels)
+                vcf_writer.write_genotypes(chromosome, variant_table, only_snvs)
                 logger.info("Done writing VCF")
 
             logger.debug("Chromosome %r finished", chromosome)
@@ -411,8 +412,9 @@ def add_arguments(parser):
         help='Reduce coverage to at most MAXCOV (default: %(default)s).')
     arg('--mapping-quality', '--mapq', metavar='QUAL',
         default=20, type=int, help='Minimum mapping quality (default: %(default)s)')
-    arg('--indels', dest='indels', default=False, action='store_true',
-        help='Also genotype indels (default: genotype only SNPs)')
+    arg('--indels', dest='indels_used', action='store_true',
+        help=SUPPRESS)
+    arg('--only-snvs', action='store_true', help='Genotype only SNVs')
     arg('--ignore-read-groups', default=False, action='store_true',
         help='Ignore read groups in BAM header and assume all reads come '
         'from the same sample.')
@@ -483,7 +485,10 @@ def validate(args, parser):
         parser.error("Option --use-ped-samples can only be used when PED file is provided (--ped).")
     if args.use_ped_samples and args.samples:
         parser.error("Option --use-ped-samples cannot be used together with --samples")
+    if args.indels_used:
+        logger.warning("Ignoring --indels as indel genotyping is default in WhatsHap 2.0+")
 
 
 def main(args):
+    del args.indels_used
     run_genotype(**vars(args))
diff --git a/whatshap/cli/haplotag.py b/whatshap/cli/haplotag.py
index cf41f3d9..5ce15685 100644
--- a/whatshap/cli/haplotag.py
+++ b/whatshap/cli/haplotag.py
@@ -498,7 +498,7 @@ def run_haplotag(
         timers.start("haplotag-init")
         try:
             vcf_reader = stack.enter_context(
-                VcfReader(variant_file, indels=True, phases=True, ploidy=ploidy)
+                VcfReader(variant_file, only_snvs=False, phases=True, ploidy=ploidy)
             )
         except OSError as err:
             raise CommandLineError(f"Error while loading variant file {variant_file}: {err}")
@@ -529,7 +529,7 @@ def run_haplotag(
                 None if reference is False else reference,
                 NumericSampleIds(),
                 ignore_read_groups,
-                indels=False,
+                only_snvs=False,
                 duplicates=True,
             )
         )
diff --git a/whatshap/cli/phase.py b/whatshap/cli/phase.py
index 61c523bb..7229244f 100644
--- a/whatshap/cli/phase.py
+++ b/whatshap/cli/phase.py
@@ -286,7 +286,7 @@ def run_whatshap(
     samples: Optional[Sequence[str]] = None,
     chromosomes: Optional[List[str]] = None,
     ignore_read_groups: bool = False,
-    indels: bool = True,
+    only_snvs: bool = False,
     mapping_quality: int = 20,
     read_merging: bool = False,
     read_merging_error_rate: float = 0.15,
@@ -374,7 +374,7 @@ def run_whatshap(
                 numeric_sample_ids,
                 ignore_read_groups,
                 mapq_threshold=mapping_quality,
-                indels=indels,
+                only_snvs=only_snvs,
             )
         )
         show_phase_vcfs = phased_input_reader.has_vcfs
@@ -393,7 +393,7 @@ def run_whatshap(
                     in_path=variant_file,
                     out_file=output,
                     tag=tag,
-                    indels=indels,
+                    only_snvs=only_snvs,
                 )
             )
         except (OSError, VcfError) as e:
@@ -401,7 +401,7 @@ def run_whatshap(
 
         # Only read genotype likelihoods from VCFs when distrusting genotypes
         vcf_reader = stack.enter_context(
-            VcfReader(variant_file, indels=indels, genotype_likelihoods=distrust_genotypes)
+            VcfReader(variant_file, only_snvs=only_snvs, genotype_likelihoods=distrust_genotypes)
         )
 
         if ignore_read_groups and not samples and len(vcf_reader.samples) > 1:
@@ -1041,8 +1041,8 @@ def add_arguments(parser):
         "quality marginally. Avoid using this in the normal case! (default: %(default)s)")
     arg("--mapping-quality", "--mapq", metavar="QUAL",
         default=20, type=int, help="Minimum mapping quality (default: %(default)s)")
-    arg("--indels", dest="indels", default=False, action="store_true",
-        help="Also phase indels (default: do not phase indels)")
+    arg("--indels", dest="indels_used", action="store_true", help=SUPPRESS)
+    arg("--only-snvs", default=False, action="store_true", help="Phase only SNVs")
     arg("--ignore-read-groups", default=False, action="store_true",
         help="Ignore read groups in BAM/CRAM header and assume all reads come "
         "from the same sample.")
@@ -1158,6 +1158,8 @@ def validate(args, parser):
             "The experimental --full-genotyping option has been removed. Instead, please run "
             "'whatshap genotype' prior to running 'whatshap phase'"
         )
+    if args.indels_used:
+        logger.warning("Ignoring --indels as indel phasing is default in WhatsHap 2.0+")
 
 
 def main(args):
@@ -1166,4 +1168,5 @@ def main(args):
     del args.no_reference
     del args.max_coverage_was_used
     del args.full_genotyping
+    del args.indels_used
     run_whatshap(**vars(args))
diff --git a/whatshap/cli/polyphase.py b/whatshap/cli/polyphase.py
index 85060590..cf62e3a7 100644
--- a/whatshap/cli/polyphase.py
+++ b/whatshap/cli/polyphase.py
@@ -48,7 +48,7 @@ def run_polyphase(
     chromosomes=None,
     verify_genotypes=False,
     ignore_read_groups=False,
-    indels=False,
+    only_snvs=False,
     mapping_quality=20,
     tag="PS",
     include_haploid_sets=False,
@@ -93,7 +93,7 @@ def run_polyphase(
                 reference,
                 numeric_sample_ids,
                 ignore_read_groups,
-                indels=indels,
+                only_snvs=only_snvs,
                 mapq_threshold=mapping_quality,
             )
         )
@@ -104,14 +104,14 @@ def run_polyphase(
         else:
             command_line = None
         try:
-            vcf_writer = stack.enter_context(
+            vcf_writer: PhasedVcfWriter = stack.enter_context(
                 PhasedVcfWriter(
                     command_line=command_line,
                     in_path=variant_file,
                     out_file=output,
                     tag=tag,
                     ploidy=ploidy,
-                    indels=indels,
+                    only_snvs=only_snvs,
                     include_haploid_sets=include_haploid_sets,
                     mav=mav,
                 )
@@ -122,7 +122,7 @@ def run_polyphase(
         vcf_reader = stack.enter_context(
             VcfReader(
                 variant_file,
-                indels=indels,
+                only_snvs=only_snvs,
                 phases=True,
                 genotype_likelihoods=False,
                 ploidy=ploidy,
@@ -417,13 +417,8 @@ def add_arguments(parser):
         type=int,
         help="Minimum mapping quality (default: %(default)s)",
     )
-    arg(
-        "--indels",
-        dest="indels",
-        default=False,
-        action="store_true",
-        help="Also phase indels (default: do not phase indels)",
-    )
+    arg("--indels", dest="indels_used", action="store_true", help=argparse.SUPPRESS)
+    arg("--only-snvs", action="store_true", help="Only phase SNVs")
     arg(
         "--ignore-read-groups",
         default=False,
@@ -545,7 +540,10 @@ def add_arguments(parser):
 def validate(args, parser):
     if args.block_cut_sensitivity > 5 or args.block_cut_sensitivity < 0:
         parser.error("Block cut sensitivity must be an integer value between 0 and 5.")
+    if args.indels_used:
+        logger.warning("Ignoring --indels as indel phasing is default in WhatsHap 2.0+")
 
 
 def main(args):
+    del args.indels_used
     run_polyphase(**vars(args))
diff --git a/whatshap/cli/polyphasegenetic.py b/whatshap/cli/polyphasegenetic.py
index c88bc62d..18b5f76e 100644
--- a/whatshap/cli/polyphasegenetic.py
+++ b/whatshap/cli/polyphasegenetic.py
@@ -65,7 +65,7 @@ def run_polyphasegenetic(
     output=sys.stdout,
     samples=None,
     chromosomes=None,
-    indels=False,
+    only_snvs=False,
     tag="PS",
     write_command_line_header=True,
     plot=False,
@@ -111,7 +111,7 @@ def run_polyphasegenetic(
         parent_reader = stack.enter_context(
             VcfReader(
                 variant_file,
-                indels=indels,
+                only_snvs=only_snvs,
                 phases=True,
                 genotype_likelihoods=False,
                 ploidy=ploidy,
@@ -124,7 +124,7 @@ def run_polyphasegenetic(
             progeny_reader = stack.enter_context(
                 VcfReader(
                     progeny_file,
-                    indels=indels,
+                    only_snvs=only_snvs,
                     phases=True,
                     genotype_likelihoods=False,
                     ploidy=ploidy,
@@ -157,7 +157,7 @@ def run_polyphasegenetic(
             complexity_support=complexity_support,
             ratio_cutoff=ratio_cutoff,
             distrust_genotypes=distrust_genotypes,
-            allow_deletions=indels,
+            allow_deletions=not only_snvs,
             plot=plot,
             output=output,
         )
@@ -537,13 +537,8 @@ def add_arguments(parser):
     )
 
     arg = parser.add_argument_group("Input pre-processing, selection, and filtering").add_argument
-    arg(
-        "--indels",
-        dest="indels",
-        default=False,
-        action="store_true",
-        help="Also phase indels (default: do not phase indels)",
-    )
+    arg("--indels", dest="indels_used", action="store_true", help=argparse.SUPPRESS)
+    arg("--only-snvs", action="store_true", help="Phase only SNVs")
     arg(
         "--sample",
         dest="samples",
@@ -644,6 +639,8 @@ def validate(args, parser):
         parser.error("Odd ploidies are not supported.")
     if args.ploidy < 2:
         parser.error("Ploidy must be at least 2.")
+    if args.indels_used:
+        logger.warning("Ignoring --indels as indel phasing is default in WhatsHap 2.0+")
 
 
 def main(args):
diff --git a/whatshap/cli/stats.py b/whatshap/cli/stats.py
index 4451b778..f7f626d1 100644
--- a/whatshap/cli/stats.py
+++ b/whatshap/cli/stats.py
@@ -486,7 +486,7 @@ def run_stats(
             gtf_file = stack.enter_context(open(gtf, "wt"))
             gtfwriter = GtfWriter(gtf_file)
 
-        vcf_reader = VcfReader(vcf, phases=True, indels=not only_snvs)
+        vcf_reader = VcfReader(vcf, phases=True, only_snvs=only_snvs)
         if len(vcf_reader.samples) == 0:
             logger.error("Input VCF does not contain any sample")
             return 1
diff --git a/whatshap/polyphase/plots.py b/whatshap/polyphase/plots.py
index 2339b315..5fa54e86 100644
--- a/whatshap/polyphase/plots.py
+++ b/whatshap/polyphase/plots.py
@@ -1174,7 +1174,7 @@ def create_genetic_plots(
         ]
         ground_truth_reader = VcfReader(
             ground_truth_file,
-            indels=False,
+            only_snvs=False,
             phases=True,
             genotype_likelihoods=False,
             ploidy=phasing_param.ploidy,
diff --git a/whatshap/vcf.py b/whatshap/vcf.py
index 8c5b7c03..df2d191d 100644
--- a/whatshap/vcf.py
+++ b/whatshap/vcf.py
@@ -496,7 +496,7 @@ class VcfReader:
     def __init__(
         self,
         path: Union[str, PathLike],
-        indels: bool = False,
+        only_snvs: bool = False,
         phases: bool = False,
         genotype_likelihoods: bool = False,
         ignore_genotypes: bool = False,
@@ -506,14 +506,13 @@ class VcfReader:
     ):
         """
         path -- Path to VCF file
-        indels -- Whether to include also insertions and deletions in the list of
-            variants.
+        only_snvs -- Whether to only include SNVs in the list of variants.
         ignore_genotypes -- In case of genotyping algorithm, no genotypes may be given in
                                 vcf, so ignore all genotypes
         ploidy -- Ploidy of the samples
         """
         # TODO Always include deletions since they can 'overlap' other variants
-        self._indels = indels
+        self._only_snvs = only_snvs
         self._vcf_reader = VariantFile(os.fspath(path))
         self._path = path
         self._phases = phases
@@ -653,7 +652,7 @@ class VcfReader:
                 n_snvs += 1
             else:
                 n_other += 1
-                if not self._indels:
+                if self._only_snvs:
                     continue
 
             if (prev_position is not None) and (prev_position > pos):
@@ -1077,7 +1076,7 @@ class PhasedVcfWriter(VcfAugmenter):
         tag: str = "PS",
         ploidy: int = 2,
         include_haploid_sets: bool = False,
-        indels: bool = False,
+        only_snvs: bool = False,
         mav: bool = False,
     ):
         """
@@ -1095,7 +1094,7 @@ class PhasedVcfWriter(VcfAugmenter):
         super().__init__(in_path, command_line, out_file, include_haploid_sets)
         self._phase_tag_found_warned = False
         self._set_phasing_tags = self._set_HP if tag == "HP" else self._set_PS
-        self._indels = indels
+        self._only_snvs = only_snvs
         self._mav = mav
 
     def setup_header(self, header: VariantHeader):
@@ -1199,8 +1198,8 @@ class PhasedVcfWriter(VcfAugmenter):
             if pos == prev_pos:
                 # duplicate position, skip it
                 continue
-            is_indel = len(str(record.ref)) > 1 or len(str(record.alts[0])) > 1
-            if not self._indels and is_indel:
+            is_snv = len(str(record.ref)) == 1 and len(str(record.alts[0])) == 1
+            if self._only_snvs and not is_snv:
                 continue
 
             # Determine whether the variant is phased in any sample
@@ -1322,7 +1321,7 @@ class GenotypeVcfWriter(VcfAugmenter):
         )
 
     def write_genotypes(
-        self, chromosome: str, variant_table: VariantTable, indels, ploidy: int = 2
+        self, chromosome: str, variant_table: VariantTable, only_snvs, ploidy: int = 2
     ) -> None:
         """
         Add genotyping information to all variants on a single chromosome.


Test patch:
diff --git a/tests/test_offspringscoring.py b/tests/test_offspringscoring.py
index cc876e29..6421c2bc 100644
--- a/tests/test_offspringscoring.py
+++ b/tests/test_offspringscoring.py
@@ -50,7 +50,7 @@ def test_correct_variant_types():
     table = list(
         VcfReader(
             "tests/data/polyphasegenetic.test.parents.vcf",
-            indels=True,
+            only_snvs=False,
             genotype_likelihoods=False,
             ploidy=4,
             mav=True,
@@ -59,7 +59,7 @@ def test_correct_variant_types():
     ptable = list(
         VcfReader(
             "tests/data/polyphasegenetic.test.progeny.vcf.gz",
-            indels=True,
+            only_snvs=False,
             genotype_likelihoods=False,
             ploidy=4,
             mav=True,
@@ -89,7 +89,7 @@ def test_compute_gt_likelihoods():
     table = list(
         VcfReader(
             "tests/data/polyphasegenetic.test.parents.vcf",
-            indels=True,
+            only_snvs=False,
             genotype_likelihoods=False,
             ploidy=4,
             mav=True,
@@ -98,7 +98,7 @@ def test_compute_gt_likelihoods():
     ptable = list(
         VcfReader(
             "tests/data/polyphasegenetic.test.progeny.vcf.gz",
-            indels=True,
+            only_snvs=False,
             genotype_likelihoods=False,
             ploidy=4,
             mav=True,
diff --git a/tests/test_polyphasegenetic.py b/tests/test_polyphasegenetic.py
index 0ce08395..b80d6b6c 100644
--- a/tests/test_polyphasegenetic.py
+++ b/tests/test_polyphasegenetic.py
@@ -6,7 +6,7 @@ def test_read_parent_vcf():
     tables = list(
         VcfReader(
             "tests/data/polyphasegenetic.test.parents.vcf",
-            indels=True,
+            only_snvs=False,
             genotype_likelihoods=False,
             ploidy=4,
             mav=True,
@@ -31,7 +31,7 @@ def test_read_progeny_vcf():
     tables = list(
         VcfReader(
             "tests/data/polyphasegenetic.test.progeny.vcf.gz",
-            indels=True,
+            only_snvs=False,
             genotype_likelihoods=False,
             ploidy=4,
             mav=True,
diff --git a/tests/test_readsetreader.py b/tests/test_readsetreader.py
index 1a7e2c62..2b83775a 100644
--- a/tests/test_readsetreader.py
+++ b/tests/test_readsetreader.py
@@ -101,10 +101,10 @@ def test_allele_dection_01():
         reference=None,
         numeric_sample_ids=NumericSampleIds(),
         ignore_read_groups=True,
-        indels=True,
+        only_snvs=False,
         mapq_threshold=20,
     )
-    vcf_reader = VcfReader(f"{path}.vcf", phases=False, indels=True)
+    vcf_reader = VcfReader(f"{path}.vcf", phases=False, only_snvs=False)
     sample = vcf_reader.samples[0]
     table = list(vcf_reader)[0]
     chromosome = table.chromosome
@@ -127,10 +127,10 @@ def test_allele_dection_02():
         reference=None,
         numeric_sample_ids=NumericSampleIds(),
         ignore_read_groups=True,
-        indels=True,
+        only_snvs=False,
         mapq_threshold=20,
     )
-    vcf_reader = VcfReader(f"{path}.vcf", phases=False, indels=True)
+    vcf_reader = VcfReader(f"{path}.vcf", phases=False, only_snvs=False)
     sample = vcf_reader.samples[0]
     table = list(vcf_reader)[0]
     chromosome = table.chromosome
@@ -154,10 +154,10 @@ def test_allele_dection_03():
         reference=None,
         numeric_sample_ids=NumericSampleIds(),
         ignore_read_groups=True,
-        indels=True,
+        only_snvs=False,
         mapq_threshold=20,
     )
-    vcf_reader = VcfReader(f"{path}.vcf", phases=False, indels=True)
+    vcf_reader = VcfReader(f"{path}.vcf", phases=False, only_snvs=False)
     sample = vcf_reader.samples[0]
     table = list(vcf_reader)[0]
     chromosome = table.chromosome
@@ -184,10 +184,10 @@ def test_allele_dection_04():
         reference=None,
         numeric_sample_ids=NumericSampleIds(),
         ignore_read_groups=True,
-        indels=True,
+        only_snvs=False,
         mapq_threshold=20,
     )
-    vcf_reader = VcfReader(f"{path}.vcf", phases=False, indels=True)
+    vcf_reader = VcfReader(f"{path}.vcf", phases=False, only_snvs=False)
     sample = vcf_reader.samples[0]
     table = list(vcf_reader)[0]
     chromosome = table.chromosome
@@ -212,10 +212,10 @@ def test_allele_dection_05():
             reference=ref,
             numeric_sample_ids=NumericSampleIds(),
             ignore_read_groups=True,
-            indels=True,
+            only_snvs=False,
             mapq_threshold=20,
         )
-        vcf_reader = VcfReader(f"{path}.vcf", phases=False, indels=True)
+        vcf_reader = VcfReader(f"{path}.vcf", phases=False, only_snvs=False)
         sample = vcf_reader.samples[0]
         table = list(vcf_reader)[0]
         chromosome = table.chromosome
@@ -242,10 +242,10 @@ def test_allele_dection_multi_01():
         reference=None,
         numeric_sample_ids=NumericSampleIds(),
         ignore_read_groups=True,
-        indels=True,
+        only_snvs=False,
         mapq_threshold=20,
     )
-    vcf_reader = VcfReader(f"{path}.vcf", phases=False, indels=True, mav=True)
+    vcf_reader = VcfReader(f"{path}.vcf", phases=False, only_snvs=False, mav=True)
     sample = vcf_reader.samples[0]
     table = list(vcf_reader)[0]
     chromosome = table.chromosome
@@ -266,10 +266,10 @@ def test_allele_dection_multi_02():
         reference="tests/data/alleledetection.fasta",
         numeric_sample_ids=NumericSampleIds(),
         ignore_read_groups=True,
-        indels=True,
+        only_snvs=False,
         mapq_threshold=20,
     )
-    vcf_reader = VcfReader(f"{path}.vcf", phases=False, indels=True, mav=True)
+    vcf_reader = VcfReader(f"{path}.vcf", phases=False, only_snvs=False, mav=True)
     sample = vcf_reader.samples[0]
     table = list(vcf_reader)[0]
     chromosome = table.chromosome
diff --git a/tests/test_run_genotype.py b/tests/test_run_genotype.py
index 7ec0ba97..c98ac538 100644
--- a/tests/test_run_genotype.py
+++ b/tests/test_run_genotype.py
@@ -80,7 +80,7 @@ def test_with_reference():
 
 
 @pytest.mark.parametrize("priors", [True, False])
-def test_no_indels(tmpdir, priors):
+def test_only_snvs(tmpdir, priors):
     prioroutput = str(tmpdir.join("priors.vcf")) if priors else None
     outvcf = str(tmpdir.join("output_gl.vcf"))
     run_genotype(
@@ -88,7 +88,7 @@ def test_no_indels(tmpdir, priors):
         variant_file="tests/data/pacbio/variants.vcf",
         reference="tests/data/pacbio/reference.fasta",
         output=outvcf,
-        indels=False,
+        only_snvs=True,
         nopriors=not priors,
         prioroutput=prioroutput,
     )
@@ -118,7 +118,7 @@ def test_multiallelic(tmpdir):
         variant_file="tests/data/multiallelic.vcf",
         reference="tests/data/pacbio/reference.fasta",
         output=outvcf,
-        indels=False,
+        only_snvs=True,
     )
     vcf_reader = VariantFile(outvcf)
     for record in vcf_reader:
@@ -152,7 +152,7 @@ def test_gt_quality_threshold(threshold, tmpdir):
         variant_file="tests/data/trio.vcf",
         output=out_vcf,
         gt_qual_threshold=threshold,
-        indels=False,
+        only_snvs=True,
         prioroutput=priors_vcf,
     )
 
@@ -445,7 +445,7 @@ def test_adding_constant(constant, tmpdir):
         variant_file="tests/data/trio.vcf",
         prioroutput=priors_raw_vcf,
         output=outvcf_raw_vcf,
-        indels=False,
+        only_snvs=True,
     )
 
     # run genotyping with modified priors
@@ -454,7 +454,7 @@ def test_adding_constant(constant, tmpdir):
         variant_file="tests/data/trio.vcf",
         prioroutput=priors_const_vcf,
         output=outvcf_const_vcf,
-        indels=False,
+        only_snvs=True,
         constant=constant,
     )
 
diff --git a/tests/test_run_phase.py b/tests/test_run_phase.py
index 2993201c..566a4327 100644
--- a/tests/test_run_phase.py
+++ b/tests/test_run_phase.py
@@ -172,7 +172,7 @@ def test_with_reference_and_indels(algorithm):
         phase_input_files=["tests/data/pacbio/pacbio.bam"],
         variant_file="tests/data/pacbio/variants.vcf",
         reference="tests/data/pacbio/reference.fasta",
-        indels=True,
+        only_snvs=False,
         algorithm=algorithm,
     )
 
@@ -558,11 +558,11 @@ def test_genetic_phasing_symbolic_alt(tmpdir):
         variant_file="tests/data/trio-symbolic-alt.vcf",
         output=outvcf,
         ped="tests/data/trio.ped",
-        indels=True,
+        only_snvs=False,
     )
     assert os.path.isfile(outvcf)
 
-    tables = list(VcfReader(outvcf, phases=True, indels=True))
+    tables = list(VcfReader(outvcf, phases=True, only_snvs=False))
     assert len(tables) == 1
     table = tables[0]
     assert table.chromosome == "1"
@@ -924,7 +924,7 @@ def test_indel_phasing(algorithm, tmp_path):
     outvcf = tmp_path / "output.vcf"
     run_whatshap(
         phase_input_files=[indels_bamfile],
-        indels=True,
+        only_snvs=False,
         variant_file="tests/data/indels.vcf",
         reference="tests/data/random0.fasta",
         output=outvcf,
@@ -932,7 +932,7 @@ def test_indel_phasing(algorithm, tmp_path):
     )
     assert os.path.isfile(outvcf)
 
-    tables = list(VcfReader(outvcf, indels=True, phases=True))
+    tables = list(VcfReader(outvcf, only_snvs=False, phases=True))
     assert len(tables) == 1
     table = tables[0]
     assert table.chromosome == "random0"
@@ -969,13 +969,13 @@ def test_vcf_with_missing_headers(algorithm):
 def test_distrust_genotypes_assertion(tmp_path):
     outvcf = tmp_path / "output.vcf"
     run_whatshap(
-        indels=False,
+        only_snvs=True,
         phase_input_files=[dist_geno_bamfile],
         variant_file="tests/data/test_dist_geno.vcf",
         output=outvcf,
     )
     assert os.path.isfile(outvcf)
-    tables = list(VcfReader(outvcf, phases=True))
+    tables = list(VcfReader(outvcf, phases=True, only_snvs=True))
     assert len(tables) == 1
     table = tables[0]
     assert table.chromosome == "chr1"
diff --git a/tests/test_run_polyphase.py b/tests/test_run_polyphase.py
index 15235b0b..d812c4bf 100644
--- a/tests/test_run_polyphase.py
+++ b/tests/test_run_polyphase.py
@@ -204,11 +204,12 @@ def test_polyphase_indels(tmp_path):
         phase_input_files=["tests/data/polyploid.indels.bam"],
         variant_file="tests/data/polyploid.indels.vcf",
         ploidy=4,
+        only_snvs=True,
         ignore_read_groups=True,
         output=outvcf,
     )
     assert os.path.isfile(outvcf)
-    tables = list(VcfReader(outvcf, phases=True, indels=True))
+    tables = list(VcfReader(outvcf, phases=True, only_snvs=False))
     table = tables[0]
 
     run_polyphase(
@@ -217,10 +218,10 @@ def test_polyphase_indels(tmp_path):
         ploidy=4,
         ignore_read_groups=True,
         output=outvcf_indel,
-        indels=True,
+        only_snvs=False,
     )
     assert os.path.isfile(outvcf_indel)
-    tables_indel = list(VcfReader(outvcf_indel, phases=True, indels=True))
+    tables_indel = list(VcfReader(outvcf_indel, phases=True, only_snvs=False))
     table_indel = tables_indel[0]
 
     sample = "Test_Indel"
diff --git a/tests/test_variantselection.py b/tests/test_variantselection.py
index 73e910cc..02de4c7d 100644
--- a/tests/test_variantselection.py
+++ b/tests/test_variantselection.py
@@ -72,7 +72,7 @@ def test_compute_phasable_variants_1():
     tables = list(
         VcfReader(
             "tests/data/polyphasegenetic.test.parents.vcf",
-            indels=True,
+            only_snvs=False,
             genotype_likelihoods=False,
             ploidy=4,
             mav=True,
@@ -141,7 +141,7 @@ def test_compute_phasable_variants_2():
     tables = list(
         VcfReader(
             "tests/data/polyphasegenetic.test.parents.vcf",
-            indels=True,
+            only_snvs=False,
             genotype_likelihoods=False,
             ploidy=4,
             mav=True,
@@ -209,7 +209,7 @@ def test_compute_phasable_variants_3():
     tables = list(
         VcfReader(
             "tests/data/polyphasegenetic.test.parents.vcf",
-            indels=True,
+            only_snvs=False,
             genotype_likelihoods=False,
             ploidy=4,
             mav=True,
diff --git a/tests/test_vcf.py b/tests/test_vcf.py
index 3e145394..74acee50 100644
--- a/tests/test_vcf.py
+++ b/tests/test_vcf.py
@@ -328,7 +328,7 @@ def test_read_duplicate_position():
     # As soon as we can actually work with multiple such rows, this test
     # needs to be updated since it currently just checks whether the second of
     # the positions is skipped.
-    table = list(VcfReader("tests/data/duplicate-positions.vcf", indels=True))[0]
+    table = list(VcfReader("tests/data/duplicate-positions.vcf", only_snvs=False))[0]
     assert len(table.variants) == 2
     assert table.variants[0].position == 1
     assert table.variants[0].reference_allele == "A"
@@ -420,7 +420,7 @@ def test_read_region():
 
 def test_read_region_subsets():
     regions = [(1069570, 1070690), (1074910, 1076152)]
-    vcf_reader = VcfReader("tests/data/haplotag_1.vcf.gz", indels=True)
+    vcf_reader = VcfReader("tests/data/haplotag_1.vcf.gz", only_snvs=False)
     table = vcf_reader.fetch_regions("chr1", regions)
     assert table.chromosome == "chr1"
     assert len(table.variants) == 8