Repository: iossifovlab/gpf

Issue:
Extract region splitting functionality into a separate util
Region splitting functionality can be found in multiple places:
- genomic_scores
- reference_genome
- annotate_vcf
- import tools

We should extract all of this functionality into a set of utilities that cover all our cases.

NOTES:
- Try extracting get_chromosome_length from GenomicPositionTable as a separate util for getting lengths of contigs not found in the reference genome, and then utilise it in our region splitting functions
- We expect to have two functions: one that creates a dictionary of (contig)->(length of contig) and another that converts such a dictionary to a set of regions, split by a given region size
- We should utilise the Region class in region_utils 

Discussion:

Patch:
diff --git a/dae/dae/genomic_resources/genomic_position_table/table.py b/dae/dae/genomic_resources/genomic_position_table/table.py
index e53617d65..340d136fd 100644
--- a/dae/dae/genomic_resources/genomic_position_table/table.py
+++ b/dae/dae/genomic_resources/genomic_position_table/table.py
@@ -165,37 +165,3 @@ class GenomicPositionTable(abc.ABC):
         This is to be overwritten by the subclass. It should return a list of
         the chromomes in the file in the order determinted by the file.
         """
-
-    def get_chromosome_length(self, chrom, step=100_000_000):
-        """Return the length of a chromosome (or contig).
-
-        Returned value is guarnteed to be larget than the actual contig length.
-        """
-        def any_records(riter):
-            try:
-                next(riter)
-            except StopIteration:
-                return False
-
-            return True
-
-        # First we find any region that includes the last record i.e.
-        # the length of the chromosome
-        left, right = None, None
-        pos = step
-        while left is None or right is None:
-            if any_records(self.get_records_in_region(chrom, pos, None)):
-                left = pos
-                pos = pos * 2
-            else:
-                right = pos
-                pos = pos // 2
-        # Second we use binary search to narrow the region until we find the
-        # index of the last element (in left) and the length (in right)
-        while (right - left) > 5_000_000:
-            pos = (left + right) // 2
-            if any_records(self.get_records_in_region(chrom, pos, None)):
-                left = pos
-            else:
-                right = pos
-        return right
diff --git a/dae/dae/genomic_resources/genomic_scores.py b/dae/dae/genomic_resources/genomic_scores.py
index 0762dbf07..e1eaad185 100644
--- a/dae/dae/genomic_resources/genomic_scores.py
+++ b/dae/dae/genomic_resources/genomic_scores.py
@@ -15,6 +15,7 @@ from functools import lru_cache
 from jinja2 import Template
 
 from dae.task_graph.graph import TaskGraph
+from dae.utils.regions import split_into_regions, get_chromosome_length_tabix
 from . import GenomicResource
 from .reference_genome import build_reference_genome_from_resource
 from .resource_implementation import GenomicResourceImplementation, \
@@ -24,6 +25,7 @@ from .resource_implementation import GenomicResourceImplementation, \
     ResourceConfigValidationMixin
 from .genomic_position_table import build_genomic_position_table, Line, \
     TabixGenomicPositionTable, VCFGenomicPositionTable, VCFLine
+from .genomic_position_table.table_inmemory import InmemoryGenomicPositionTable
 from .histogram import Histogram, HistogramStatisticMixin
 from .statistics.min_max import MinMaxValue, MinMaxValueStatisticMixin
 
@@ -646,12 +648,21 @@ class GenomicScore(
         return ref_genome
 
     def _get_chrom_regions(self, region_size, grr=None):
+        regions = []
         ref_genome_id = self.resource.get_labels().get("reference_genome")
         ref_genome = self._get_reference_genome_cached(grr, ref_genome_id)
-        if ref_genome is not None:
-            regions = self._split_into_regions(region_size, ref_genome)
-        else:
-            regions = self._split_into_regions(region_size)
+        for chrom in self.get_all_chromosomes():
+            if ref_genome is not None:
+                chrom_length = ref_genome.get_chrom_length(chrom)
+            else:
+                if isinstance(self.table, InmemoryGenomicPositionTable):
+                    raise ValueError("In memory tables are not supported")
+                chrom_length = get_chromosome_length_tabix(
+                    self.table.pysam_file, chrom
+                )
+            regions.extend(
+                split_into_regions(chrom, chrom_length, region_size)
+            )
         return regions
 
     def _add_min_max_tasks(self, graph, region_size, grr=None):
@@ -663,7 +674,10 @@ class GenomicScore(
         """
         min_max_tasks = []
         regions = self._get_chrom_regions(region_size, grr)
-        for chrom, start, end in regions:
+        for region in regions:
+            chrom = region.chrom
+            start = region.start
+            end = region.stop
             min_max_tasks.append(graph.create_task(
                 f"{self.score_id}_calculate_min_max_{chrom}_{start}_{end}",
                 GenomicScore._do_min_max,
@@ -740,7 +754,10 @@ class GenomicScore(
         """
         histogram_tasks = []
         regions = self._get_chrom_regions(region_size, grr)
-        for chrom, start, end in regions:
+        for region in regions:
+            chrom = region.chrom
+            start = region.start
+            end = region.stop
             histogram_tasks.append(graph.create_task(
                 f"{self.score_id}_calculate_histogram_{chrom}_{start}_{end}",
                 GenomicScore._do_histogram,
@@ -835,34 +852,6 @@ class GenomicScore(
                 score_histogram.plot(outfile)
         return merged_histograms
 
-    def _split_into_regions(self, region_size, reference_genome=None):
-        chromosomes = self.get_all_chromosomes()
-        for chrom in chromosomes:
-            if reference_genome is not None \
-                    and chrom in reference_genome.chromosomes:
-                chrom_len = reference_genome.get_chrom_length(chrom)
-            else:
-                if reference_genome is not None:
-                    logger.info(
-                        "chromosome %s of %s not found in reference genome %s",
-                        chrom, self.resource.resource_id,
-                        reference_genome.resource_id
-                    )
-                else:
-                    logger.info(
-                        "chromosome %s of %s using table, no reference genome",
-                        chrom, self.resource.resource_id
-                    )
-                chrom_len = self.table.get_chromosome_length(chrom)
-            logger.debug(
-                "Chromosome '%s' has length %s",
-                chrom, chrom_len)
-            i = 1
-            while i < chrom_len - region_size:
-                yield chrom, i, i + region_size - 1
-                i += region_size
-            yield chrom, i, None
-
     def calc_info_hash(self):
         """Compute and return the info hash."""
         return "infohash"
diff --git a/dae/dae/genomic_resources/reference_genome.py b/dae/dae/genomic_resources/reference_genome.py
index fc573be71..814a3328b 100644
--- a/dae/dae/genomic_resources/reference_genome.py
+++ b/dae/dae/genomic_resources/reference_genome.py
@@ -471,7 +471,14 @@ class ReferenceGenome(
             yield chrom, i, None
 
     def fetch(self, chrom, start, stop, buffer_size=512):
-        """Yield the nucleotides in a specific region."""
+        """
+        Yield the nucleotides in a specific region.
+
+        While line feed calculation can be inaccurate because not every fetch
+        will start at the start of a line, line feeds add extra characters
+        to read and the output is limited by the amount of nucleotides
+        expected to be read.
+        """
         if chrom not in self.chromosomes:
             logger.warning(
                 "chromosome %s not found in %s",
@@ -496,41 +503,19 @@ class ReferenceGenome(
         total_length = length + line_feeds
         read_progress = 0
 
-        while read_progress < total_length:
+        while read_progress < length:
             read_length = min(buffer_size, total_length - read_progress)
-            line_feeds = 1 + read_length // self._index[chrom]["seqLineLength"]
             sequence = self._sequence.read(read_length).decode("ascii")
             sequence = sequence.replace("\n", "").upper()
-            sequence = sequence[:read_length - line_feeds]
+            end = min(read_progress + read_length, length - read_progress)
+            sequence = sequence[:end]
             for nuc in sequence:
                 yield nuc
-            read_progress += read_length
-
-    def get_sequence_new(self, chrom, start, stop):
-        """Return sequence of nucleotides from specified chromosome region."""
-        return "".join(self.fetch(chrom, start, stop))
+            read_progress += len(sequence)
 
     def get_sequence(self, chrom, start, stop):
         """Return sequence of nucleotides from specified chromosome region."""
-        if chrom not in self.chromosomes:
-            logger.warning(
-                "chromosome %s not found in %s",
-                chrom, self.resource.resource_id)
-            return None
-        assert self._sequence is not None
-        self._sequence.seek(
-            self._index[chrom]["startBit"]
-            + start
-            - 1
-            + (start - 1) // self._index[chrom]["seqLineLength"]
-        )
-
-        length = stop - start + 1
-        line_feeds = 1 + length // self._index[chrom]["seqLineLength"]
-
-        sequence = self._sequence.read(length + line_feeds).decode("ascii")
-        sequence = sequence.replace("\n", "")[:length]
-        return sequence.upper()
+        return "".join(self.fetch(chrom, start, stop))
 
     def pair_iter(self, chrom, start, stop):
         """
diff --git a/dae/dae/utils/regions.py b/dae/dae/utils/regions.py
index 0b62e8a10..948cf1b7f 100644
--- a/dae/dae/utils/regions.py
+++ b/dae/dae/utils/regions.py
@@ -2,6 +2,8 @@ from __future__ import annotations
 from collections import defaultdict
 import copy
 
+import pysam
+
 
 def bedfile2regions(bed_filename):
     """Transform BED file into list of regions."""
@@ -25,6 +27,62 @@ def regions2bedfile(regions, bed_filename):
                 f"{reg.chrom}\t{reg.start-1}\t{reg.stop}\n")
 
 
+def split_into_regions(
+        chrom: str, chrom_length: int,
+        region_size: int) -> list[Region]:
+    """Return a list of regions for a chrom with a given length."""
+    regions = []
+
+    current_start = 1
+    while current_start < chrom_length + 1:
+        end = min(chrom_length, current_start + region_size - 1)
+        regions.append(Region(chrom, current_start, end))
+        current_start = current_start + region_size
+
+    return regions
+
+
+def get_chromosome_length_tabix(
+    tabix_file: pysam.TabixFile, chrom: str,
+    step=100_000_000, precision=5_000_000
+):
+    """
+    Return the length of a chromosome (or contig).
+
+    Returned value is guarnteed to be larget than the actual contig length.
+    """
+    def any_records(riter):
+        try:
+            next(riter)
+        except StopIteration:
+            return False
+
+        return True
+
+    # First we find any region that includes the last record i.e.
+    # the length of the chromosome
+    left, right = None, None
+    pos = step
+    while left is None or right is None:
+        region = Region(chrom, pos, None)
+        if any_records(tabix_file.fetch(str(region))):
+            left = pos
+            pos = pos * 2
+        else:
+            right = pos
+            pos = pos // 2
+    # Second we use binary search to narrow the region until we find the
+    # index of the last element (in left) and the length (in right)
+    while (right - left) > precision:
+        pos = (left + right) // 2
+        region = Region(chrom, pos, None)
+        if any_records(tabix_file.fetch(str(region))):
+            left = pos
+        else:
+            right = pos
+    return right
+
+
 class Region:
     """Class representing a genomic region."""
 


Test patch:
diff --git a/dae/dae/genomic_resources/tests/test_cli_stats.py b/dae/dae/genomic_resources/tests/test_cli_stats.py
index 85d9dc96d..a2d595a34 100644
--- a/dae/dae/genomic_resources/tests/test_cli_stats.py
+++ b/dae/dae/genomic_resources/tests/test_cli_stats.py
@@ -9,7 +9,7 @@ from dae.genomic_resources.repository import GR_CONF_FILE_NAME
 from dae.genomic_resources.histogram import Histogram
 from dae.genomic_resources.testing import \
     setup_directories, convert_to_tab_separated, \
-    build_filesystem_test_repository
+    build_filesystem_test_repository, setup_tabix
 from dae.genomic_resources.resource_implementation import \
     GenomicResourceImplementation, ResourceStatistics
 
@@ -117,7 +117,8 @@ def test_stats_allele_score(tmp_path):
             GR_CONF_FILE_NAME: """
                 type: allele_score
                 table:
-                    filename: data.mem
+                    filename: data.txt.gz
+                    format: tabix
                 scores:
                     - id: freq
                       type: float
@@ -130,21 +131,23 @@ def test_stats_allele_score(tmp_path):
                       max: 1
                       x_scale: linear
                       y_scale: log
-                """,
-            "data.mem": convert_to_tab_separated("""
-                chrom  pos_begin  reference  alternative  freq
-                1      10         A          G            0.02
-                1      10         A          C            0.03
-                1      10         A          A            0.04
-                1      16         CA         G            0.03
-                1      16         C          T            0.04
-                1      16         C          A            0.05
-                2      16         CA         G            0.03
-                2      16         C          T            EMPTY
-                2      16         C          A            0.05
-            """)
+                """
         }
     })
+    setup_tabix(
+        tmp_path / "one" / "data.txt.gz",
+        """
+        #chrom  pos_begin  reference  alternative  freq
+        1      10         A          G            0.02
+        1      10         A          C            0.03
+        1      10         A          A            0.04
+        1      16         CA         G            0.03
+        1      16         C          T            0.04
+        1      16         C          A            0.05
+        2      16         CA         G            0.03
+        2      16         C          T            EMPTY
+        2      16         C          A            0.05
+        """, seq_col=0, start_col=1, end_col=1)
 
     repo = build_filesystem_test_repository(tmp_path)
 
@@ -182,7 +185,8 @@ def test_stats_position_score(tmp_path):
             GR_CONF_FILE_NAME: """
                 type: position_score
                 table:
-                    filename: data.mem
+                    filename: data.txt.gz
+                    format: tabix
                 scores:
                     - id: phastCons100way
                       type: float
@@ -209,17 +213,19 @@ def test_stats_position_score(tmp_path):
                       max: 4
                       x_scale: linear
                       y_scale: linear
-                """,
-            "data.mem": convert_to_tab_separated("""
-                chrom  pos_begin  pos_end  s1    s2
-                1      10         15       0.02  -1
-                1      17         19       0.03  0
-                1      22         25       0.46  EMPTY
-                2      5          80       0.01  3
-                2      10         11       0.02  3
-                """)
-        },
+                """
+        }
     })
+    setup_tabix(
+        tmp_path / "one" / "data.txt.gz",
+        """
+        #chrom  pos_begin  pos_end  s1    s2
+        1      10         15       0.02  -1
+        1      17         19       0.03  0
+        1      22         25       0.46  EMPTY
+        2      5          80       0.01  3
+        2      10         11       0.02  3
+        """, seq_col=0, start_col=1, end_col=2)
 
     repo = build_filesystem_test_repository(tmp_path)
 
@@ -277,7 +283,8 @@ def test_stats_np_score(tmp_path):
             GR_CONF_FILE_NAME: """
                 type: np_score
                 table:
-                    filename: data.mem
+                    filename: data.txt.gz
+                    format: tabix
                 scores:
                     - id: cadd_raw
                       type: float
@@ -303,21 +310,23 @@ def test_stats_np_score(tmp_path):
                       max: 4
                       x_scale: linear
                       y_scale: linear
-            """,
-            "data.mem": convert_to_tab_separated("""
-                chrom  pos_begin  pos_end  reference  alternative  s1    s2
-                1      10         15       A          G            0.02  2
-                1      10         15       A          C            0.03  -1
-                1      10         15       A          T            0.04  4
-                1      16         19       C          G            0.03  3
-                1      16         19       C          T            0.04  EMPTY
-                1      16         19       C          A            0.05  0
-                2      16         19       C          A            0.03  3
-                2      16         19       C          T            0.04  3
-                2      16         19       C          G            0.05  4
-            """)
+            """
         }
     })
+    setup_tabix(
+        tmp_path / "one" / "data.txt.gz",
+        """
+        #chrom  pos_begin  pos_end  reference  alternative  s1    s2
+        1      10         15       A          G            0.02  2
+        1      10         15       A          C            0.03  -1
+        1      10         15       A          T            0.04  4
+        1      16         19       C          G            0.03  3
+        1      16         19       C          T            0.04  EMPTY
+        1      16         19       C          A            0.05  0
+        2      16         19       C          A            0.03  3
+        2      16         19       C          T            0.04  3
+        2      16         19       C          G            0.05  4
+        """, seq_col=0, start_col=1, end_col=2)
 
     repo = build_filesystem_test_repository(tmp_path)
 
@@ -374,7 +383,8 @@ def test_minmax(tmp_path):
             GR_CONF_FILE_NAME: """
                 type: position_score
                 table:
-                    filename: data.mem
+                    filename: data.txt.gz
+                    format: tabix
                 scores:
                     - id: phastCons100way
                       type: float
@@ -386,19 +396,21 @@ def test_minmax(tmp_path):
                       bins: 100
                       x_scale: linear
                       y_scale: linear
-            """,
-            "data.mem": convert_to_tab_separated("""
-                chrom  pos_begin  pos_end  s1
-                1      10         15       0.0
-                1      17         19       0.03
-                1      22         25       0.46
-                2      5          80       0.01
-                2      10         11       1.0
-                3      5          17       1.0
-                3      18         20       0.01
-            """)
+            """
         }
     })
+    setup_tabix(
+        tmp_path / "one" / "data.txt.gz",
+        """
+        #chrom  pos_begin  pos_end  s1
+        1      10         15       0.0
+        1      17         19       0.03
+        1      22         25       0.46
+        2      5          80       0.01
+        2      10         11       1.0
+        3      5          17       1.0
+        3      18         20       0.01
+        """, seq_col=0, start_col=1, end_col=2)
 
     repo = build_filesystem_test_repository(tmp_path)
 
@@ -424,7 +436,8 @@ def test_reference_genome_usage(tmp_path, mocker):
             GR_CONF_FILE_NAME: """
                 type: position_score
                 table:
-                    filename: data.mem
+                    filename: data.txt.gz
+                    format: tabix
                 scores:
                     - id: phastCons100way
                       type: float
@@ -439,17 +452,7 @@ def test_reference_genome_usage(tmp_path, mocker):
                 meta:
                     labels:
                         reference_genome: genome
-            """,
-            "data.mem": convert_to_tab_separated("""
-                chrom  pos_begin  pos_end  s1
-                1      10         15       0.0
-                1      17         19       0.03
-                1      22         25       0.46
-                2      5          8        0.01
-                2      10         11       1.0
-                3      5          17       1.0
-                3      18         20       0.01
-            """)
+            """
         },
         "genome": {
             GR_CONF_FILE_NAME: """
@@ -478,6 +481,18 @@ def test_reference_genome_usage(tmp_path, mocker):
 
         }
     })
+    setup_tabix(
+        tmp_path / "one" / "data.txt.gz",
+        """
+        #chrom  pos_begin  pos_end  s1
+        1      10         15       0.0
+        1      17         19       0.03
+        1      22         25       0.46
+        2      5          8        0.01
+        2      10         11       1.0
+        3      5          17       1.0
+        3      18         20       0.01
+        """, seq_col=0, start_col=1, end_col=2)
 
     repo = build_filesystem_test_repository(tmp_path)
 
@@ -504,8 +519,7 @@ def test_reference_genome_usage(tmp_path, mocker):
 
     genomic_table_length_mock = mocker.Mock(return_value=30)
     mocker.patch(
-        "dae.genomic_resources.genomic_position_table.table"
-        ".GenomicPositionTable.get_chromosome_length",
+        "dae.genomic_resources.genomic_scores.get_chromosome_length_tabix",
         new=genomic_table_length_mock
     )
 
diff --git a/dae/dae/genomic_resources/tests/test_genomic_position_table.py b/dae/dae/genomic_resources/tests/test_genomic_position_table.py
index 7ca6c1bdd..361094ea8 100644
--- a/dae/dae/genomic_resources/tests/test_genomic_position_table.py
+++ b/dae/dae/genomic_resources/tests/test_genomic_position_table.py
@@ -11,6 +11,8 @@ from dae.genomic_resources.genomic_position_table import \
     VCFGenomicPositionTable, \
     build_genomic_position_table
 
+from dae.utils.regions import get_chromosome_length_tabix
+
 from dae.genomic_resources.testing import \
     build_inmemory_test_resource, build_filesystem_test_resource, \
     setup_directories, convert_to_tab_separated, setup_tabix, setup_vcf
@@ -1070,6 +1072,7 @@ def test_tabix_max_buffer(
         )
 
 
+@pytest.mark.xfail(reason="Unsupported method")
 def test_contig_length():
     res = build_inmemory_test_resource({
         "genomic_resource.yaml": """
@@ -1092,7 +1095,7 @@ def test_contig_length():
 
 
 def test_contig_length_tabix_table(tabix_table):
-    assert tabix_table.get_chromosome_length("1") >= 13
+    assert get_chromosome_length_tabix(tabix_table.pysam_file, "1") >= 13
 
 
 def test_vcf_autodetect_format(vcf_res_autodetect_format):
diff --git a/dae/dae/genomic_resources/tests/test_reference_genome_resource.py b/dae/dae/genomic_resources/tests/test_reference_genome_resource.py
index e62fa90a4..f7187eed5 100644
--- a/dae/dae/genomic_resources/tests/test_reference_genome_resource.py
+++ b/dae/dae/genomic_resources/tests/test_reference_genome_resource.py
@@ -210,6 +210,44 @@ def test_reference_genome_fetch(genome_fixture):
         ]
 
 
+def test_reference_genome_fetch_corner_case(genome_fixture):
+    res = build_filesystem_test_resource(genome_fixture)
+    reference_genome = build_reference_genome_from_resource(res)
+
+    with reference_genome.open():
+        result = list(reference_genome.fetch("pesho", 1, 9))
+        assert result == [
+            "N",
+            "N",
+            "A",
+            "C",
+            "C",
+            "C",
+            "A",
+            "A",
+            "A",
+        ]
+
+
+def test_reference_genome_fetch_small_buffer(genome_fixture):
+    res = build_filesystem_test_resource(genome_fixture)
+    reference_genome = build_reference_genome_from_resource(res)
+
+    with reference_genome.open():
+        result = list(reference_genome.fetch("pesho", 1, 9, 1))
+        assert result == [
+            "N",
+            "N",
+            "A",
+            "C",
+            "C",
+            "C",
+            "A",
+            "A",
+            "A",
+        ]
+
+
 def test_reference_genome_pair_iter(genome_fixture):
     res = build_filesystem_test_resource(genome_fixture)
     reference_genome = build_reference_genome_from_resource(res)
diff --git a/dae/dae/utils/tests/test_region.py b/dae/dae/utils/tests/test_region.py
index 4cec8e6ae..96b4ae33c 100644
--- a/dae/dae/utils/tests/test_region.py
+++ b/dae/dae/utils/tests/test_region.py
@@ -1,5 +1,10 @@
+# pylint: disable=W0621,C0114,C0116,W0212,W0613
 import pytest
-from dae.utils.regions import Region, collapse, collapse_no_chrom
+import pysam
+from dae.utils.regions import Region, collapse, collapse_no_chrom, \
+    split_into_regions, get_chromosome_length_tabix
+
+from dae.genomic_resources.testing import setup_tabix
 
 
 @pytest.mark.parametrize(
@@ -56,3 +61,69 @@ def test_collapse_no_chrom_simple(regions, expected):
     assert len(result) == len(expected)
     for res, exp in zip(result, expected):
         assert res == exp
+
+
+@pytest.mark.parametrize(
+    "chrom,chrom_length,region_size,expected",
+    [
+        ("1", 10, 5, [Region("1", 1, 5), Region("1", 6, 10)]),
+        ("1", 10, 100, [Region("1", 1, 10)]),
+        ("1", 4, 1, [
+            Region("1", 1, 1),
+            Region("1", 2, 2),
+            Region("1", 3, 3),
+            Region("1", 4, 4)
+        ]),
+    ]
+)
+def test_split_into_regions(chrom, chrom_length, region_size, expected):
+    result = split_into_regions(chrom, chrom_length, region_size)
+
+    assert len(result) == len(expected)
+
+    for i, res in enumerate(result):
+        assert res == expected[i], f"{res} != {expected[i]}"
+
+
+@pytest.fixture()
+def sample_tabix(tmp_path):
+    filepath = tmp_path / "data.txt.gz"
+    setup_tabix(
+        filepath,
+        """
+            #chrom pos_begin pos_end
+            1     10        12
+            1     11        11
+            1     12        13
+            1     13        14
+            2     1         2
+            2     1         4
+            3     500000    500010
+        """, seq_col=0, start_col=1, end_col=2
+    )
+    return pysam.TabixFile(str(filepath))
+
+
+@pytest.mark.parametrize(
+    "precision",
+    [
+        5_000_000,
+        1000,
+        100,
+        1,
+    ]
+)
+def test_get_chrom_length(sample_tabix, precision):
+    one_length = get_chromosome_length_tabix(
+        sample_tabix, "1", precision=precision)
+    assert one_length > 14
+    assert one_length - 14 <= precision
+    two_length = get_chromosome_length_tabix(
+        sample_tabix, "2", precision=precision)
+    assert two_length > 4
+    assert two_length - 14 <= precision
+    three_length = get_chromosome_length_tabix(
+        sample_tabix, "3", precision=precision
+    )
+    assert three_length > 500_010
+    assert three_length - 500_010 <= precision