Skip to content

Commit 5f0d4c0

Browse files
authored
Fix #190: Genome.merged_gene_intervals returns union of gene loci (#347)
Promotes the existing private _merge_ranges helper in transcript.py to a public common.merge_intervals (the old name is kept as an alias). Adds Genome.merged_gene_intervals(contig, strand=None) which queries all (start, end) pairs from the gene table on the given contig and returns a sorted list of non-overlapping inclusive intervals. Adjacent intervals (end+1 == next start) are merged. Useful for "is this position inside any gene?" / gene-density without double-counting overlapping genes. Bump version to 2.9.4.
1 parent 92e4297 commit 5f0d4c0

5 files changed

Lines changed: 135 additions & 19 deletions

File tree

pyensembl/common.py

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,25 @@ def _memoize_cache_key(args, kwargs):
5050
return tuple(cache_key_list)
5151

5252

53+
def merge_intervals(ranges):
54+
"""
55+
Sort ``[(start, end)]`` inclusive-inclusive ranges and merge any pair
56+
that overlaps or is exactly adjacent (``end + 1 == next start``).
57+
Returns a new list of merged ``(start, end)`` tuples in ascending order.
58+
"""
59+
if not ranges:
60+
return []
61+
ordered = sorted(ranges)
62+
merged = [ordered[0]]
63+
for start, end in ordered[1:]:
64+
prev_start, prev_end = merged[-1]
65+
if start <= prev_end + 1:
66+
merged[-1] = (prev_start, max(prev_end, end))
67+
else:
68+
merged.append((start, end))
69+
return merged
70+
71+
5372
def memoize(fn):
5473
"""Simple reset-able memoization decorator for functions and methods,
5574
assumes that all arguments to the function can be hashed and

pyensembl/genome.py

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,9 +22,11 @@
2222
from serializable import Serializable
2323

2424
from .download_cache import DownloadCache
25+
from .common import merge_intervals
2526
from .database import Database
2627
from .exon import Exon
2728
from .gene import Gene
29+
from .normalization import normalize_chromosome, normalize_strand
2830
from .search import find_nearest_locus
2931
from .sequence_data import SequenceData
3032
from .transcript import Transcript
@@ -588,6 +590,24 @@ def nearest_gene(self, contig, position, end=None, strand=None):
588590
loci=self.genes(contig=contig, strand=strand),
589591
)
590592

593+
def merged_gene_intervals(self, contig, strand=None):
594+
"""
595+
Return the union of all gene loci on ``contig`` as a sorted list of
596+
non-overlapping ``(start, end)`` tuples. Adjacent intervals
597+
(``end + 1 == next start``) are merged into one. Pass ``strand`` to
598+
restrict to one strand of the contig.
599+
600+
Useful for asking "is this position inside any gene?" or computing
601+
gene-density coverage without double-counting overlapping genes.
602+
"""
603+
sql = "SELECT start, end FROM gene WHERE seqname = ?"
604+
params = [normalize_chromosome(contig)]
605+
if strand is not None:
606+
sql += " AND strand = ?"
607+
params.append(normalize_strand(strand))
608+
rows = self.db.run_sql_query(sql, query_params=params)
609+
return merge_intervals(rows)
610+
591611
def nearest_transcript(self, contig, position, end=None, strand=None):
592612
"""
593613
Find the transcript on ``contig`` whose locus is nearest to the

pyensembl/transcript.py

Lines changed: 4 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -12,27 +12,13 @@
1212

1313
from memoized_property import memoized_property
1414

15-
from .common import memoize
15+
from .common import memoize, merge_intervals
1616
from .exon import Exon
1717
from .locus_with_genome import LocusWithGenome
1818

1919

20-
def _merge_ranges(ranges):
21-
"""
22-
Sort [(start, end)] inclusive-inclusive ranges and merge any that are
23-
adjacent or overlapping (end+1 == next start).
24-
"""
25-
if not ranges:
26-
return []
27-
ordered = sorted(ranges)
28-
merged = [ordered[0]]
29-
for start, end in ordered[1:]:
30-
prev_start, prev_end = merged[-1]
31-
if start <= prev_end + 1:
32-
merged[-1] = (prev_start, max(prev_end, end))
33-
else:
34-
merged.append((start, end))
35-
return merged
20+
# Back-compat alias for callers that imported the private helper.
21+
_merge_ranges = merge_intervals
3622

3723

3824
class Transcript(LocusWithGenome):
@@ -435,7 +421,7 @@ def coding_sequence_position_ranges(self):
435421
ranges.extend(
436422
self._transcript_feature_position_ranges("stop_codon", required=False)
437423
)
438-
return _merge_ranges(ranges)
424+
return merge_intervals(ranges)
439425

440426
@memoized_property
441427
def complete(self):

pyensembl/version.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
__version__ = "2.9.3"
1+
__version__ = "2.9.4"
22

33
def print_version():
44
print(f"v{__version__}")
Lines changed: 91 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,91 @@
1+
"""
2+
Issue #190: Genome.merged_gene_intervals(contig) returns sorted
3+
non-overlapping (start, end) tuples after merging adjacent/overlapping
4+
gene loci on the contig.
5+
"""
6+
7+
from pyensembl.common import merge_intervals
8+
9+
from .common import TemporaryDirectory, eq_
10+
from pyensembl import Genome
11+
12+
13+
def test_merge_intervals_handles_overlap_adjacency_and_disjoint():
14+
# Overlapping
15+
eq_(merge_intervals([(1, 10), (5, 15)]), [(1, 15)])
16+
# Adjacent (end+1 == next start) merges
17+
eq_(merge_intervals([(1, 10), (11, 20)]), [(1, 20)])
18+
# One-base gap stays separate
19+
eq_(merge_intervals([(1, 10), (12, 20)]), [(1, 10), (12, 20)])
20+
# Unsorted input is reordered
21+
eq_(merge_intervals([(50, 60), (1, 10), (5, 15)]), [(1, 15), (50, 60)])
22+
# Containment
23+
eq_(merge_intervals([(1, 100), (10, 20)]), [(1, 100)])
24+
# Empty
25+
eq_(merge_intervals([]), [])
26+
27+
28+
def test_merged_gene_intervals_against_synthetic_gtf():
29+
"""Build a tiny GTF with three overlapping and one disjoint gene on
30+
chr1 plus a single gene on chr2, then assert the merged intervals."""
31+
with TemporaryDirectory() as tmpdir:
32+
import os
33+
gtf_path = os.path.join(tmpdir, "merge.gtf")
34+
with open(gtf_path, "w") as f:
35+
# chr1: gene1 [100,200], gene2 [150,250] overlap → [100,250]
36+
# gene3 [251,300] adjacent → merges into [100,300]
37+
# gene4 [500,600] disjoint
38+
# chr1 -: gene5 on minus strand for strand-filter test
39+
# chr2: gene6 standalone
40+
for gene_id, contig, start, end, strand in [
41+
("G1", "1", 100, 200, "+"),
42+
("G2", "1", 150, 250, "+"),
43+
("G3", "1", 251, 300, "+"),
44+
("G4", "1", 500, 600, "+"),
45+
("G5", "1", 700, 800, "-"),
46+
("G6", "2", 1000, 1100, "+"),
47+
]:
48+
f.write(
49+
"%s\ttest\tgene\t%d\t%d\t.\t%s\t.\t"
50+
'gene_id "%s"; gene_name "%s";\n'
51+
% (contig, start, end, strand, gene_id, gene_id)
52+
)
53+
f.write(
54+
"%s\ttest\ttranscript\t%d\t%d\t.\t%s\t.\t"
55+
'gene_id "%s"; transcript_id "%s_t"; gene_name "%s";\n'
56+
% (contig, start, end, strand, gene_id, gene_id, gene_id)
57+
)
58+
f.write(
59+
"%s\ttest\texon\t%d\t%d\t.\t%s\t.\t"
60+
'gene_id "%s"; transcript_id "%s_t"; exon_number "1"; gene_name "%s";\n'
61+
% (contig, start, end, strand, gene_id, gene_id, gene_id)
62+
)
63+
64+
genome = Genome(
65+
reference_name="GRCh38",
66+
annotation_name="merge_intervals_test",
67+
gtf_path_or_url=gtf_path,
68+
cache_directory_path=tmpdir,
69+
)
70+
genome.index()
71+
72+
# chr1 (any strand): overlapping triplet collapses, gene4 stays,
73+
# minus-strand gene5 is included when strand is not filtered.
74+
eq_(
75+
genome.merged_gene_intervals("1"),
76+
[(100, 300), (500, 600), (700, 800)],
77+
)
78+
# Strand-restricted: gene5 is excluded
79+
eq_(
80+
genome.merged_gene_intervals("1", strand="+"),
81+
[(100, 300), (500, 600)],
82+
)
83+
# Only the minus-strand gene
84+
eq_(
85+
genome.merged_gene_intervals("1", strand="-"),
86+
[(700, 800)],
87+
)
88+
# chr2: single gene
89+
eq_(genome.merged_gene_intervals("2"), [(1000, 1100)])
90+
# Empty contig returns []
91+
eq_(genome.merged_gene_intervals("3"), [])

0 commit comments

Comments
 (0)