More work on filtering. Added MAPQS to format line. MAPQ can be filte…

…red using dysgu filter. Option to set search depth when SV calling
kcleal · Mar 22, 2024 · 3faafbe · 3faafbe
1 parent 952d43f
commit 3faafbe
Show file tree

Hide file tree

Showing 7 changed files with 123 additions and 84 deletions.
diff --git a/dysgu/call_component.pyx b/dysgu/call_component.pyx
@@ -772,19 +772,19 @@ cdef single(rds, int insert_size, int insert_stdev, float insert_ppf, int clip_l
             svlen_adjusted = int(np.median([b[0] for b in size_pos_bounded]))
             posA_adjusted = int(np.median([b[1] for b in size_pos_bounded]))
             posB_adjusted = int(np.median([b[2] for b in size_pos_bounded]))
-            if svtype_m == 2:  # del
-                posA = posA_adjusted
-                posB = posB_adjusted
-                svlen = svlen_adjusted
-                er.preciseA = True
-                er.preciseB = True
-            else:
+            # if svtype_m == 2:  # del
+            #     posA = posA_adjusted
+            #     posB = posB_adjusted
+            #     svlen = svlen_adjusted
+            #     er.preciseA = True
+            #     er.preciseB = True
+            # else:
                 # to_assemble = False
-                er.preciseA = True
-                er.preciseB = True
-                posA = posA_adjusted
-                posB = posB_adjusted
-                svlen = svlen_adjusted
+            er.preciseA = True
+            er.preciseB = True
+            posA = posA_adjusted
+            posB = posB_adjusted
+            svlen = svlen_adjusted
 
         else:
             svlen = int(np.median([sp[4] for sp in spanning_alignments]))

diff --git a/dysgu/cluster.pyx b/dysgu/cluster.pyx
@@ -837,6 +837,7 @@ def pipe1(args, infile, kind, regions, ibam, ref_genome, sample_name, bam_iter=N
                                             temp_dir=tdir,
                                             find_n_aligned_bases=find_n_aligned_bases,
                                             position_distance_thresh=args['sd'],
+                                            max_search_depth=args['search_depth']
     )
     sites_index = None
     if sites_adder:

diff --git a/dysgu/extra_metrics.pyx b/dysgu/extra_metrics.pyx
@@ -262,7 +262,7 @@ cpdef filter_poorly_aligned_ends(spanning_alignments, float divergence=0.02):
     return spanning, 1 - (len(spanning) / len(spanning_alignments))
 
 
-cpdef gap_size_upper_bound(AlignedSegment alignment, int cigarindex, int pos_input, int end_input, int length_extend=15, float divergence=0.02):
+cpdef gap_size_upper_bound(AlignedSegment alignment, int cigarindex, int pos_input, int end_input, int length_extend=15, float divergence=0.02, float len_t=3):
     # expand indel using cigar, merges nearby gaps into the SV event
     cdef int pos, end, l, opp, extent_left, extent_right, candidate_type, candidate_len, len_input, i, dist, dist_thresh, middle, last_seen_size
     cdef uint32_t cigar_value
@@ -276,7 +276,7 @@ cpdef gap_size_upper_bound(AlignedSegment alignment, int cigarindex, int pos_inp
     candidate_type = cigar_value & BAM_CIGAR_MASK
     candidate_len = cigar_value >> BAM_CIGAR_SHIFT
     len_input = candidate_len
-    dist_thresh = min(len_input * 3, <int> (<float> log2_32(1 + candidate_len) / divergence))
+    dist_thresh = min(<int> (len_input * len_t), <int> (<float> log2_32(1 + candidate_len) / divergence))
     i = cigarindex + 1
     dist = 0
     last_seen_size = candidate_len
@@ -308,7 +308,7 @@ cpdef gap_size_upper_bound(AlignedSegment alignment, int cigarindex, int pos_inp
                 else:
                     candidate_len += l
                     extent_right = end
-            dist_thresh = min(l * 3, <int> (<float> log2_32(1 + candidate_len) / divergence))
+            dist_thresh = min(<int> (l * len_t), <int> (<float> log2_32(1 + candidate_len) / divergence))
             dist = 0
         i += 1
 
@@ -342,7 +342,7 @@ cpdef gap_size_upper_bound(AlignedSegment alignment, int cigarindex, int pos_inp
                 else:
                     candidate_len += l
                     extent_left = pos
-            dist_thresh = min(l * 3, <int>( <float>log2_32(1 + candidate_len) / divergence) )
+            dist_thresh = min(<int> (l * len_t), <int>( <float>log2_32(1 + candidate_len) / divergence) )
             dist = 0
         i -= 1