modify bugs for singe breakend functions

friend1ws · Mar 8, 2021 · 5eac94b · 5eac94b
1 parent 020ed8f
commit 5eac94b
Show file tree

Hide file tree

Showing 3 changed files with 13 additions and 3 deletions.
diff --git a/nanomonsv/cluster_sbnd.py b/nanomonsv/cluster_sbnd.py
@@ -139,7 +139,7 @@ def cluster_supporting_reads_sbnd(input_file, output_file, control_bed = None,
         for line in hin:
 
             tchr, _, tpos, treadid, _, tdir, tinfo = line.rstrip('\n').split('\t')
-            tstart = int(tpos) - sbnd_cluster_margin_size
+            tstart = max(int(tpos) - sbnd_cluster_margin_size, 0)
             tend = int(tpos) + sbnd_cluster_margin_size
 
             # skip single breakend from secondary alignment

diff --git a/nanomonsv/locate_bp.py b/nanomonsv/locate_bp.py
@@ -17,9 +17,15 @@ def get_refined_bp(contig, fasta_file_ins, chr1, start1, end1, dir1, chr2, start
 
     if mode != "i":
 
+        ref_len1 = fasta_file_ins.get_reference_length(chr1)
+        ref_len2 = fasta_file_ins.get_reference_length(chr2)
+
         bstart1, bend1 = max(int(start1) - rd_margin, 1), int(end1) + rd_margin
         bstart2, bend2 = max(int(start2) - rd_margin, 1), int(end2) + rd_margin
 
+        if ref_len1 < bend1: bend1 = ref_len1
+        if ref_len2 < bend2: bend2 = ref_len2    
+
         region1_seq = fasta_file_ins.fetch(chr1, bstart1 - 1, bend1)
         region2_seq = fasta_file_ins.fetch(chr2, bstart2 - 1, bend2)
 

diff --git a/nanomonsv/locate_bp_sbnd.py b/nanomonsv/locate_bp_sbnd.py
@@ -15,10 +15,14 @@ def get_refined_bp_sbnd(tconsensus, fasta_file_h, tchr, tstart, tend, tdir, hout
 
     tconsensus_part = tconsensus[:1000] if len(tconsensus) > 1000 else tconsensus
 
+    ref_len = fasta_file_h.get_reference_length(tchr)
+    if tstart < 1: tstart = 1
+    if ref_len < tend: tend = ref_len 
+
     if tdir == '+':
-        qseq = fasta_file_h.fetch(tchr, int(tstart) - margin, int(tend))
+        qseq = fasta_file_h.fetch(tchr, max(int(tstart) - margin, 0), int(tend))
     else:
-        qseq = fasta_file_h.fetch(tchr, int(tstart), int(tend) + margin)
+        qseq = fasta_file_h.fetch(tchr, max(int(tstart) - 1, 0), int(tend) + margin)
         qseq = reverse_complement(qseq)
 
     user_matrix = parasail.matrix_create("ACGT", 1, -2)