From 63699ea95946b31c7732c5bfbb737ca53c0530b5 Mon Sep 17 00:00:00 2001
From: Bertache <samir.bertache@ens-lyon.fr>
Date: Tue, 25 Mar 2025 17:05:21 +0100
Subject: [PATCH] Doctest valide et complet
---
hi-classifier/classify.py | 3 +-
hi-classifier/filter_binary.py | 207 ++++++++++++++++++++++++++-------
hi-classifier/positionfrag.py | 198 ++++++++++++++++++++++++++++++-
hi-classifier/readmanager.py | 22 ++--
4 files changed, 374 insertions(+), 56 deletions(-)
diff --git a/hi-classifier/classify.py b/hi-classifier/classify.py
index 73fded3..f69d72b 100644
--- a/hi-classifier/classify.py
+++ b/hi-classifier/classify.py
@@ -9,10 +9,9 @@ from bam import read_bam_pair
from digest import build_restriction_index_from_records
from processing import process_items
from processmanager import ProcessManager
+from restrictionsite import SearchInDataBase
from writer import write_bam_pair
-from restriction_site import SearchInDataBase
-
def communicate(write_processes, compute_processes):
print(
diff --git a/hi-classifier/filter_binary.py b/hi-classifier/filter_binary.py
index 1d9e32b..15b06b8 100644
--- a/hi-classifier/filter_binary.py
+++ b/hi-classifier/filter_binary.py
@@ -303,13 +303,13 @@ def condition_dangling_left(paired_read, TOLERANCE):
>>> read1 = ("r1", 0, "chr1", 105, 30, "50M", None, None, 0, "ACGT")
>>> read2 = ("r2", 0, "chr1", 200, 30, "50M", None, None, 0, "ACGT")
>>> pr = PairedRead(read1, read2)
- >>> pr.re_proximal_1 = 100 # site de restriction à 100
+ >>> pr.re_proximal_1 = (1, 100) # site de restriction à 100
>>> condition_dangling_left(pr, TOLERANCE=10)
True
>>> condition_dangling_left(pr, TOLERANCE=4)
False
"""
- return abs(paired_read.start_1 - paired_read.re_proximal_1) <= TOLERANCE
+ return abs(paired_read.start_1 - paired_read.re_proximal_1[1]) <= TOLERANCE
def condition_dangling_right(paired_read, TOLERANCE):
@@ -321,13 +321,13 @@ def condition_dangling_right(paired_read, TOLERANCE):
>>> read1 = ("r1", 0, "chr1", 100, 30, "50M", None, None, 0, "ACGT")
>>> read2 = ("r2", 0, "chr1", 208, 30, "50M", None, None, 0, "ACGT")
>>> pr = PairedRead(read1, read2)
- >>> pr.re_proximal_2 = 210
+ >>> pr.re_proximal_2 = (1, 210)
>>> condition_dangling_right(pr, TOLERANCE=5)
True
>>> condition_dangling_right(pr, TOLERANCE=1)
False
"""
- return abs(paired_read.start_2 - paired_read.re_proximal_2) <= TOLERANCE
+ return abs(paired_read.start_2 - paired_read.re_proximal_2[1]) <= TOLERANCE
def condition_dangling_both(paired_read, TOLERANCE):
@@ -336,11 +336,11 @@ def condition_dangling_both(paired_read, TOLERANCE):
Examples:
>>> from readmanager import PairedRead
- >>> read1 = ("r1", 0, "chr1", 103, 30, "50M", None, None, 0, "ACGT")
- >>> read2 = ("r2", 0, "chr1", 107, 30, "50M", None, None, 0, "ACGT")
+ >>> read1 = ("r1", 0, "chr1", 105, 30, "50M", None, None, 0, "ACGT")
+ >>> read2 = ("r2", 16, "chr1", 208, 30, "50M", None, None, 0, "ACGT")
>>> pr = PairedRead(read1, read2)
- >>> pr.re_proximal_1 = 100
- >>> pr.re_proximal_2 = 105
+ >>> pr.re_proximal_1 = (1, 100)
+ >>> pr.re_proximal_2 = (2, 210)
>>> condition_dangling_both(pr, TOLERANCE=5)
True
>>> condition_dangling_both(pr, TOLERANCE=1)
@@ -541,20 +541,145 @@ def Classified(paired_read, tolerance=3):
Examples:
>>> from readmanager import PairedRead
- >>> # Imaginons un simple "Contact valide" type "Close"
- >>> read1 = ("r1", 0, "chr1", 100, 30, "50M", None, None, 0, "ACGT")
- >>> read2 = ("r2", 0x10, "chr1", 300, 30, "50M", None, None, 0, "ACGT") # brin négatif
+ >>> from positionfrag import annotate_intra_paired_read
+ >>> restr_index = {
+ ... "chr1": {
+ ... "sites": [
+ ... (1, 100),
+ ... (2, 200),
+ ... (3, 300),
+ ... (4, 400)
+ ... ],
+ ... "fragments": [
+ ... (1, (0, 100)),
+ ... (2, (100, 200)),
+ ... (3, (200, 300)),
+ ... (4, (300, 400)),
+ ... (5, (400, 500))
+ ... ]
+ ... }
+ ... }
+ >>> read1 = ("r1", 0, "chr1", 101, 30, "50M", None, None, 0, "ACGT")
+ >>> read2 = ("r2", 0, "chr1", 300, 30, "50M", None, None, 0, "ACGT") # brin négatif
+ >>> pr = PairedRead(read1, read2)
+ >>> # On annote pr avec restr_index
+ >>> annotate_intra_paired_read(pr, restr_index)
+ >>> pr.re_1
+ (2, 200)
+ >>> pr.strand_1
+ '+'
+ >>> pr.start_1
+ 101
+ >>> pr.end_1
+ 151
+ >>> pr.fragment_1
+ [(2, (100, 200))]
+ >>> pr.re_2
+ (4, 400)
+ >>> pr.strand_2
+ '+'
+ >>> pr.start_2
+ 300
+ >>> pr.end_2
+ 350
+ >>> pr.fragment_2
+ [(4, (300, 400))]
+ >>> # On appelle maintenant Classified
+ >>> label = Classified(pr, tolerance=5)
+ >>> label
+ 'Up'
+ >>> #################################################################
+ >>> # Exemple 1 : "Close"
+ >>> # Conditions : re_distant True, orientation_divergent True,
+ >>> # single_fragment True, frag_identical False, frag_adjacent False.
+ >>> read1 = ("r1", 16, "chr1", 256, 30, "50M", None, None, 0, "ACGT")
+ >>> read2 = ("r2", 0, "chr1", 326, 30, "50M", None, None, 0, "ACGT") # read2 sur '-'
>>> pr = PairedRead(read1, read2)
- >>> pr.fragment_1 = [(10, (100,150))]
- >>> pr.fragment_2 = [(12, (300,350))]
- >>> pr.re_1 = (5, 105)
- >>> pr.re_2 = (10, 310)
- >>> # => condition_distant_RE doit être True (5 vs 10 => ecart=5>1)
- >>> # orientation_divergent => True
- >>> # single_fragment => True => ni identique ni adjacent => (10 vs 12 => ecart=2)
+ >>> annotate_intra_paired_read(pr, restr_index)
>>> label = Classified(pr, tolerance=5)
>>> label
'Close'
+ >>> #################################################################
+ >>> # Exemple 2 : "Up"
+ >>> read1 = ("r1", 0, "chr1", 42, 30, "50M", None, None, 0, "ACGT")
+ >>> read2 = ("r2", 0, "chr1", 265, 30, "50M", None, None, 0, "ACGT")
+ >>> pr = PairedRead(read1, read2)
+ >>> annotate_intra_paired_read(pr, restr_index)
+ >>> label = Classified(pr, tolerance=5)
+ >>> label
+ 'Up'
+ >>> #################################################################
+ >>> # Exemple 3 : "Far"
+ >>> read1 = ("r1", 0, "chr1", 42, 30, "50M", None, None, 0, "ACGT")
+ >>> read2 = ("r2", 18, "chr1", 358, 30, "50M", None, None, 0, "ACGT")
+ >>> pr = PairedRead(read1, read2)
+ >>> annotate_intra_paired_read(pr, restr_index)
+ >>> label = Classified(pr, tolerance=5)
+ >>> label
+ 'Far'
+ >>> #################################################################
+ # Exemple 4 : "Joined-Dangling-Left"
+ >>> read1 = ("r1", 0, "chr1", 102, 30, "50M", None, None, 0, "ACGT")
+ >>> read2 = ("r2", 16, "chr1", 291, 30, "70M", None, None, 0, "ACGT")
+ >>> pr = PairedRead(read1, read2)
+ >>> annotate_intra_paired_read(pr, restr_index)
+ >>> label = Classified(pr, tolerance=4)
+ >>> label
+ 'Joined-Dangling-Left'
+ >>> label = Classified(pr, tolerance=1)
+ >>> label
+ 'Joined'
+ >>> label = Classified(pr, tolerance=10)
+ >>> label
+ 'Joined-Dangling-Both'
+ >>> #################################################################
+ >>> # Exemple 5 : "Self-Circle"
+ >>> read1 = ("r1", 48, "chr1", 140, 30, "32M", None, None, 0, "ACGT")
+ >>> read2 = ("r2", 0, "chr1", 151, 30, "45M", None, None, 0, "ACGT")
+ >>> pr = PairedRead(read1, read2)
+ >>> annotate_intra_paired_read(pr, restr_index)
+ >>> label = Classified(pr, tolerance=5)
+ >>> label
+ 'Self-Circle'
+ >>> #################################################################
+ # Exemple 6 : "Joined" (informative)
+ >>> read1 = ("r1", 0, "chr1", 120, 30, "50M", None, None, 0, "ACGT")
+ >>> read2 = ("r2", 16, "chr1", 254, 30, "50M", None, None, 0, "ACGT")
+ >>> pr = PairedRead(read1, read2)
+ >>> annotate_intra_paired_read(pr, restr_index)
+ >>> label = Classified(pr, tolerance=5)
+ >>> label
+ 'Joined'
+ >>> #################################################################
+ >>> # Exemple 7 : "DanglingLeft"
+ >>> read1 = ("r1", 0, "chr1", 110, 30, "50M", None, None, 0, "ACGT")
+ >>> read2 = ("r2", 16, "chr1", 185, 30, "30M", None, None, 0, "ACGT")
+ >>> pr = PairedRead(read1, read2)
+ >>> annotate_intra_paired_read(pr, restr_index)
+ >>> label = Classified(pr, tolerance=1)
+ >>> label
+ 'Other'
+ >>> label = Classified(pr, tolerance=10)
+ >>> label
+ 'DanglingLeft'
+ >>> #################################################################
+ >>> # Exemple 8 : "DanglingBoth"
+ >>> read1 = ("r1", 0, "chr1", 100, 30, "50M", None, None, 0, "ACGT")
+ >>> read2 = ("r2", 16, "chr1", 195, 30, "50M", None, None, 0, "ACGT")
+ >>> pr = PairedRead(read1, read2)
+ >>> annotate_intra_paired_read(pr, restr_index)
+ >>> label = Classified(pr, tolerance=5)
+ >>> label
+ 'DanglingBoth'
+ >>> #################################################################
+ >>> # Exemple 9 : "Other"
+ >>> read1 = ("r1", 0, "chr1", 100, 30, "50M", None, None, 0, "ACGT")
+ >>> read2 = ("r2", 0, "chr1", 100, 30, "50M", None, None, 0, "ACGT")
+ >>> pr = PairedRead(read1, read2)
+ >>> annotate_intra_paired_read(pr, restr_index)
+ >>> label = Classified(pr, tolerance=5)
+ >>> label
+ 'Other'
"""
mask = compute_mask(paired_read, tolerance)
@@ -597,8 +722,8 @@ def Classified(paired_read, tolerance=3):
return "Down"
return "Down"
- # Si on ne matche aucune orientation "typique", on reste en "ValideOther"
- return "ValideOther"
+ # Si on ne matche aucune orientation "possible", on passe en "Erreur_Algorythmique"
+ return "Erreur_Algorythmique" # Amas Zingue !
# ---------------------------------------------------------
# 2) CONTACT INFORMATIF (re_distant == False)
@@ -618,12 +743,21 @@ def Classified(paired_read, tolerance=3):
if mask["multiple_fragment"] and (
mask["frag_multiple_adjacent"] or mask["frag_multiple_identical"]
):
- return "Joined-Multiple"
+ return "Joined"
# Si single => adjacent
if mask["single_fragment"] and mask["frag_adjacent"]:
- return "Joined-Single"
- # Sinon => "Joined" par défaut
- return "Joined"
+ # Combinaison spécifique :
+ if mask["dangling_both"]:
+ return "Joined-Dangling-Both"
+
+ elif mask["dangling_left"]:
+ return "Joined-Dangling-Left"
+
+ elif mask["dangling_right"]:
+ return "Joined-Dangling-Right"
+
+ else:
+ return "Joined"
# ---------------------------------------------------------
# 3) SEMI-BIOTINYLÉS : single + frag_identical + orientation_convergent
@@ -640,31 +774,16 @@ def Classified(paired_read, tolerance=3):
db = mask["dangling_both"] # généralement = dl & dr
if db:
- return "Semi-Biot DanglingBoth"
+ return "DanglingBoth"
elif dl:
- return "Semi-Biot DanglingLeft"
+ return "DanglingLeft"
elif dr:
- return "Semi-Biot DanglingRight"
+ return "DanglingRight"
else:
- return "Semi-Biot"
-
- # ---------------------------------------------------------
- # 4) COMBINAISONS (Joined × Dangling) - exemple
- # ---------------------------------------------------------
- # Si on veut repérer un contact "Joined" qui en plus a un "dangling":
- # Dans la logique ci-dessus, "Joined" est déjà rendu plus haut.
- # On pourrait faire un second test si on veut détecter la combinaison
- # qui n'a pas été prise en compte. Ex:
- if (mask["orientation_convergent"] and mask["re_identical"]) and (
- mask["dangling_left"]
- or mask["dangling_right"]
- or mask["dangling_both"]
- ):
- # => "Joined + Dangling"
- return "Joined-Dangling"
+ return "Other" # :o
# ---------------------------------------------------------
- # 5) DIVERS => "Other"
+ # 4) DIVERS => "Other"
# ---------------------------------------------------------
return "Other"
diff --git a/hi-classifier/positionfrag.py b/hi-classifier/positionfrag.py
index 02bf16c..00e6f0f 100644
--- a/hi-classifier/positionfrag.py
+++ b/hi-classifier/positionfrag.py
@@ -2,11 +2,106 @@ from bisect import bisect_left
def get_re_site_for_read(chrom, pos, strand, restr_index):
+ """
+ Pour une position 'pos' sur le chromosome 'chrom' et pour un brin donné (strand: '+' ou '-'),
+ renvoie le tuple (site_num, position) du site de restriction le plus proche dans la direction du read,
+ en ignorant le site si sa position est exactement égale à pos (la première base du read).
+
+ - Pour un read forward ('+'): retourne le premier site dont la position est > pos.
+ - Pour un read reverse ('-'): retourne le dernier site dont la position est < pos.
+
+ Si aucun site n'est trouvé, retourne None.
+
+ Examples:
+ >>> restr_index = {
+ ... "chr1": {
+ ... "sites": [
+ ... (1, 100),
+ ... (2, 200),
+ ... (3, 300),
+ ... (4, 400),
+ ... (5, 500)
+ ... ],
+ ... "fragments": []
+ ... }
+ ... }
+ >>> get_re_site_for_read("chr1", 50, "+", restr_index)
+ (1, 100)
+ >>> get_re_site_for_read("chr1", 250, "+", restr_index)
+ (3, 300)
+ >>> get_re_site_for_read("chr1", 450, "+", restr_index)
+ (5, 500)
+ >>> get_re_site_for_read("chr1", 200, "-", restr_index)
+ (1, 100)
+ >>> get_re_site_for_read("chr1", 450, "-", restr_index)
+ (4, 400)
+ """
+ sites = restr_index[chrom]["sites"] # Liste de tuples (site_num, position)
+ if not sites:
+ return None
+ positions = [s[1] for s in sites]
+ if strand == "+":
+ idx = bisect_left(positions, pos)
+ if idx < len(positions):
+ if positions[idx] == pos:
+ # Si le site trouvé est exactement à pos, on retourne le suivant (s'il existe)
+ if idx < len(positions) - 1:
+ return sites[idx + 1]
+ else:
+ return None
+ return sites[idx]
+ else:
+ return sites[-1]
+ elif strand == "-":
+ idx = bisect_left(positions, pos)
+ if idx == 0:
+ return None
+ else:
+ if positions[idx - 1] == pos:
+ if idx - 1 > 0:
+ return sites[idx - 2]
+ else:
+ return None
+ return sites[idx - 1]
+ else:
+ raise ValueError("Orientation inconnue : " + strand)
+
+
+def get_proximal_re_site_for_read(chrom, pos, strand, restr_index):
"""
Pour une position 'pos' sur le chromosome 'chrom' et pour un brin donné (strand: '+' ou '-'),
renvoie le tuple (site_num, position) du site de restriction le plus proche dans la direction du read.
- Pour un read forward ('+'): retourne le premier site dont la position est >= pos.
- Pour un read reverse ('-'): retourne le dernier site dont la position est <= pos.
+
+ Examples:
+ >>> # On définit un index de restriction fictif
+ >>> restr_index = {
+ ... "chr1": {
+ ... "sites": [
+ ... (1, 100),
+ ... (2, 200),
+ ... (3, 300),
+ ... (4, 400)
+ ... ],
+ ... "fragments": [
+ ... (1, (0, 100)),
+ ... (2, (100, 200)),
+ ... (3, (200, 300)),
+ ... (4, (300, 400))
+ ... ]
+ ... }
+ ... }
+ >>> # Test forward ('+')
+ >>> get_re_site_for_read("chr1", 250, "+", restr_index)
+ (3, 300)
+ >>> get_re_site_for_read("chr1", 450, "+", restr_index)
+ (4, 400)
+ >>> # Test reverse ('-')
+ >>> get_re_site_for_read("chr1", 250, "-", restr_index)
+ (2, 200)
+ >>> get_re_site_for_read("chr1", 450, "-", restr_index)
+ (4, 400)
"""
sites = restr_index[chrom]["sites"] # Liste de tuples (site_num, position)
if not sites:
@@ -37,16 +132,43 @@ def get_re_proximal(chrom, pos, strand, restr_index):
AntiSens = site de restriction pour le brin "-" (reverse)
Renvoie le tuple (site_num, position) dont la position est la plus proche de pos.
+
+ Examples:
+ >>> restr_index = {
+ ... "chr1": {
+ ... "sites": [
+ ... (1, 100),
+ ... (2, 200),
+ ... (3, 300),
+ ... (4, 400)
+ ... ],
+ ... "fragments": [
+ ... (1, (0, 100)),
+ ... (2, (100, 200)),
+ ... (3, (200, 300)),
+ ... (4, (300, 400))
+ ... ]
+ ... }
+ ... }
+ >>> # pos=210 => côté forward => (3,300), distance=90
+ >>> # côté reverse => (2,200), distance=10
+ >>> # => On attend (2,200) comme site le plus proche
+ >>> get_re_proximal("chr1", 210, "+", restr_index)
+ (2, 200)
+ >>> # pos=250 => forward => (3,300), distance=50
+ >>> # reverse => (2,200), distance=50 => égalité => on renvoie AntiSens
+ >>> get_re_proximal("chr1", 240, "+", restr_index)
+ (2, 200)
"""
- Sens = get_re_site_for_read(chrom, pos, "+", restr_index)
- AntiSens = get_re_site_for_read(chrom, pos, "-", restr_index)
+ Sens = get_proximal_re_site_for_read(chrom, pos, "+", restr_index)
+ AntiSens = get_proximal_re_site_for_read(chrom, pos, "-", restr_index)
if Sens is None:
return AntiSens
if AntiSens is None:
return Sens
- if abs(Sens - pos) < abs(AntiSens - pos):
+ if abs(Sens[1] - pos) < abs(AntiSens[1] - pos):
return Sens
else:
return AntiSens
@@ -57,6 +179,32 @@ def get_fragments_for_read(chrom, start, end, restr_index):
Renvoie la liste de fragments (chaque fragment étant un tuple (frag_num, (frag_start, frag_end)))
dans lesquels se trouve la lecture.
L'intervalle de la lecture est défini par [min(start, end), max(start, end)).
+
+ Il faut plus d'une base qui se chevauche pour qu'on considère le fragment
+
+ Examples:
+ >>> restr_index = {
+ ... "chr1": {
+ ... "sites": [
+ ... (1, 100),
+ ... (2, 200),
+ ... (3, 300),
+ ... (4, 400)
+ ... ],
+ ... "fragments": [
+ ... (1, (0, 100)),
+ ... (2, (100, 200)),
+ ... (3, (200, 300)),
+ ... (4, (300, 400))
+ ... ]
+ ... }
+ ... }
+ >>> # Lecture recouvre [150, 250], donc chevauche fragments 2 (100,200) et 3 (200,300)
+ >>> get_fragments_for_read("chr1", 150, 250, restr_index)
+ [(2, (100, 200)), (3, (200, 300))]
+ >>> # Lecture recouvre [50, 99], chevauche fragment 1 (0,100)
+ >>> get_fragments_for_read("chr1", 99, 50, restr_index)
+ [(1, (0, 100))]
"""
low = min(start, end)
high = max(start, end)
@@ -77,6 +225,50 @@ def annotate_intra_paired_read(paired_read, restr_index):
le plus proche dans la direction du read.
- fragment_1 et fragment_2 : liste des fragments (peut être un ou plusieurs) couvrant l'intervalle du read.
On se base sur le chromosome, la position de début et la position de fin (déjà calculée) de chaque lecture.
+
+ NB : Intervalle à demi-ouvert tel que [start, end)
+
+ Examples:
+ >>> # On simule un index + un PairedRead
+ >>> from readmanager import PairedRead, flag_to_strand
+ >>> restr_index = {
+ ... "chr1": {
+ ... "sites": [
+ ... (1, 100),
+ ... (2, 200),
+ ... (3, 300),
+ ... (4, 400)
+ ... ],
+ ... "fragments": [
+ ... (1, (0, 100)),
+ ... (2, (100, 200)),
+ ... (3, (200, 300)),
+ ... (4, (300, 400))
+ ... ]
+ ... }
+ ... }
+ >>> read1 = ("r1", 0, "chr1", 149, 30, "50M", None, None, 0, "ACGT")
+ >>> read2 = ("r2", 16, "chr1", 351, 30, "50M", None, None, 0, "ACGT") # brin négatif
+ >>> pr = PairedRead(read1, read2)
+ >>> annotate_intra_paired_read(pr, restr_index)
+ >>> pr.re_1, pr.re_2
+ ((2, 200), (3, 300))
+ >>> pr.re_2
+ (3, 300)
+ >>> pr.re_proximal_1, pr.re_proximal_2
+ ((1, 100), (4, 400))
+ >>> pr.fragment_1
+ [(2, (100, 200))]
+ >>> pr.fragment_2
+ [(4, (300, 400))]
+ >>> read1 = ("r1", 18, "chr1", 149, 30, "50M", None, None, 0, "ACGT") # brin négatif
+ >>> read2 = ("r2", 16, "chr1", 349, 30, "50M", None, None, 0, "ACGT") # brin négatif
+ >>> pr = PairedRead(read1, read2)
+ >>> annotate_intra_paired_read(pr, restr_index)
+ >>> pr.fragment_1
+ [(1, (0, 100)), (2, (100, 200))]
+ >>> pr.fragment_2
+ [(3, (200, 300)), (4, (300, 400))]
"""
# Pour read1
chrom1 = paired_read.chrom_1
diff --git a/hi-classifier/readmanager.py b/hi-classifier/readmanager.py
index 60c946f..d4409ed 100644
--- a/hi-classifier/readmanager.py
+++ b/hi-classifier/readmanager.py
@@ -40,6 +40,10 @@ def cigar_reference_length(cigar_str: str) -> int:
0
>>> cigar_reference_length("4M")
4
+ >>> cigar_reference_length("4X15M")
+ 19
+ >>> cigar_reference_length("8H4M8H")
+ 4
"""
if not cigar_str:
return 0
@@ -68,6 +72,10 @@ def flag_to_strand(flag):
'+'
>>> flag_to_strand(20) # 20 = 0x14 => 0x10 (negatif) + 0x4 (unmapped)
'-'
+ >>> flag_to_strand(18) # 20 = 0x14 => 0x10 (negatif) + 0x4 (unmapped)
+ '-'
+ >>> flag_to_strand(12) # 20 = 0x14 => 0x10 (negatif) + 0x4 (unmapped)
+ '+'
"""
return "-" if (flag & 0x10) else "+"
@@ -104,7 +112,7 @@ class PairedRead:
>>> read2 = ("r2", 16, "chr1", 200, 30, "10M", None, None, 0, "ACGT")
>>> pr = PairedRead(read1, read2)
>>> pr
- <PairedRead r1/r2 chr=chr1 strand=(+,-) pos=(100-149, 200-191)>
+ <PairedRead r1/r2 chr=chr1 strand=(+,-) pos=(100-150, 200-190)>
>>> pr.interchrom
False
>>> # Vérifions l'ordre: read1 est plus amont => on n'a pas besoin de swap
@@ -114,7 +122,7 @@ class PairedRead:
>>> read2b = ("r2b", 16, "chr1", 100, 30, "10M", None, None, 0, "ACGT")
>>> prb = PairedRead(read1b, read2b)
>>> prb
- <PairedRead r2b/r1b chr=chr1 strand=(-,+) pos=(100-91, 500-549)>
+ <PairedRead r2b/r1b chr=chr1 strand=(-,+) pos=(100-90, 500-550)>
>>> # On voit que read2b/r1b ont été inversés pour que la coord la plus à gauche soit first.
"""
@@ -153,19 +161,19 @@ class PairedRead:
ref_len_1 = cigar_reference_length(self.cigar_1)
if self.strand_1 == "+":
# Brin positif => s'étend à droite
- self.end_1 = self.start_1 + ref_len_1 - 1
+ self.end_1 = self.start_1 + ref_len_1
else:
# Brin négatif => s'étend à gauche
- self.end_1 = self.start_1 - (ref_len_1 - 1)
+ self.end_1 = self.start_1 - ref_len_1
# --- Calcul end_2 selon self.strand_2 ---
self.end_2 = None
if self.start_2 is not None and self.cigar_2:
ref_len_2 = cigar_reference_length(self.cigar_2)
if self.strand_2 == "+":
- self.end_2 = self.start_2 + ref_len_2 - 1
+ self.end_2 = self.start_2 + ref_len_2
else:
- self.end_2 = self.start_2 - (ref_len_2 - 1)
+ self.end_2 = self.start_2 - ref_len_2
# --- Imposer l'ordre si même chromosome ---
if (
@@ -235,7 +243,7 @@ class PairedRead:
>>> read2 = ("r2", 16, "chr1", 200, 30, "50M", None, None, 0, "ACGT")
>>> pr = PairedRead(read1, read2)
>>> pr
- <PairedRead r1/r2 chr=chr1 strand=(+,-) pos=(100-149, 200-151)>
+ <PairedRead r1/r2 chr=chr1 strand=(+,-) pos=(100-150, 200-150)>
"""
return (
f"<PairedRead {self.query_name_1}/{self.query_name_2} "
--
GitLab