separated BND and INS/DEL signal scoring

Author: vinzenzmay

src/svirlpool/candidateregions/signalstrength_to_crs.py

@@ -395,14 +395,19 @@ def process_chromosome_to_proto_crs(args_tuple) -> tuple[str, Path, dict]:
     """Worker function to process one chromosome and create proto-CRs.
 
     Args:
-        args_tuple: (chr_name, chr_signals_file, chr_repeats_file, tmp_dir, buffer_region_radius)
+        args_tuple: (chr_name, chr_signals_file, chr_repeats_file, tmp_dir, buffer_region_radius, bnd_region_radius)
 
     Returns:
         (chr_name, proto_crs_file, statistics_dict)
     """
-    chr_name, chr_signals_file, chr_repeats_file, tmp_dir, buffer_region_radius = (
-        args_tuple
-    )
+    (
+        chr_name,
+        chr_signals_file,
+        chr_repeats_file,
+        tmp_dir,
+        buffer_region_radius,
+        bnd_region_radius,
+    ) = args_tuple
 
     logger.info(f"Processing chromosome {chr_name}...")
 
@@ -432,7 +437,7 @@ def process_chromosome_to_proto_crs(args_tuple) -> tuple[str, Path, dict]:
         for s in signals:
             sv_type = s[4].sv_type
             if sv_type in (3, 4):
-                margin = 3 * buffer_region_radius
+                margin = bnd_region_radius
             elif sv_type in (1, 2):
                 margin = max(
                     50, int(abs(s[4].size) * 0.5)
@@ -775,6 +780,7 @@ def create_candidate_regions(
     dropped: Path | None = None,
     bedgraph: Path | None = None,
     tmp_dir_path: Path | None = None,
+    bnd_region_radius: int = 300,
 ) -> None:
     csv.field_size_limit(sys.maxsize)
 
@@ -872,6 +878,7 @@ def create_candidate_regions(
                 chr_repeat_files.get(chr_name, tmp_dir / f"{chr_name}_repeats.bed"),
                 tmp_dir,
                 buffer_region_radius,
+                bnd_region_radius,
             ))
 
         chr_results = []
@@ -972,6 +979,7 @@ def run(args, **kwargs):
         min_cr_size=args.min_cr_size,
         dropped=args.dropped,
         tmp_dir_path=getattr(args, "tmp_dir", None),
+        bnd_region_radius=args.bnd_region_radius,
     )
 
 
@@ -1030,6 +1038,13 @@ def get_parser():
         default=300,
         help="Radius of the initial seeds' sizes around signal that can constitute a candidate region seed.",
     )
+    parser.add_argument(
+        "--bnd-region-radius",
+        type=int,
+        required=False,
+        default=300,
+        help="Margin (bp) added around BND signals when seeding candidate regions. Independent of --buffer-region-radius.",
+    )
     parser.add_argument(
         "--cutoff-median-readcount-per-region",
         type=float,

src/svirlpool/localassembly/consensus.py

@@ -518,9 +518,10 @@ def _rank_reads_by_similarity(
 
 
 def _rank_reads_by_alignment_lengths(
-        pairwise_alignment_lengths_matrix: np.ndarray,
-        sim_read_names: list[str],
-        cluster_read_names: list[str]) -> list[str]:
+    pairwise_alignment_lengths_matrix: np.ndarray,
+    sim_read_names: list[str],
+    cluster_read_names: list[str],
+) -> list[str]:
     """
     rank all reads of the cluster by the sum of their alignment lengths to all other reads in the cluster. Largest sum is ranked highest.
     """
@@ -531,7 +532,9 @@ def _rank_reads_by_alignment_lengths(
     cluster_indices = np.array([name_to_idx[r] for r in cluster_read_names])
 
     # Sub-matrix for the cluster
-    sub_sim = pairwise_alignment_lengths_matrix[np.ix_(cluster_indices, cluster_indices)]
+    sub_sim = pairwise_alignment_lengths_matrix[
+        np.ix_(cluster_indices, cluster_indices)
+    ]
     # sum alignment lengths of each read to all other reads in the cluster
     for i in range(sub_sim.shape[0]):
         sub_sim[i, i] = 0
@@ -541,6 +544,7 @@ def _rank_reads_by_alignment_lengths(
     # return the readnames sorted by sum of alignment lengths
     return [cluster_read_names[i] for i in sim_ranks]
 
+
 def find_representative_read(
     similarity_matrix: np.ndarray,
     pairwise_alignment_lengths_matrix: np.ndarray,
@@ -549,14 +553,21 @@ def find_representative_read(
 ) -> str:
     """Find the representative read of a cluster based on similarity and alignment lengths."""
     # rank reads by similarity and alignment lengths
-    ranked_by_similarity = _rank_reads_by_similarity(similarity_matrix, sim_read_names, cluster_read_names)
-    ranked_by_alignment_lengths = _rank_reads_by_alignment_lengths(pairwise_alignment_lengths_matrix, sim_read_names, cluster_read_names)
+    ranked_by_similarity = _rank_reads_by_similarity(
+        similarity_matrix, sim_read_names, cluster_read_names
+    )
+    ranked_by_alignment_lengths = _rank_reads_by_alignment_lengths(
+        pairwise_alignment_lengths_matrix, sim_read_names, cluster_read_names
+    )
     ranksums: np.ndarray = np.zeros(len(cluster_read_names))
     for i, read in enumerate(cluster_read_names):
-        ranksums[i] = ranked_by_similarity.index(read) + ranked_by_alignment_lengths.index(read)
+        ranksums[i] = ranked_by_similarity.index(
+            read
+        ) + ranked_by_alignment_lengths.index(read)
     representative_read = cluster_read_names[np.argmin(ranksums)]
     return representative_read
 
+
 def make_consensus_with_racon_subsampled(
     reference_fasta: Path,
     sam_alignments: Path,
@@ -1171,7 +1182,9 @@ def partition_reads_spectral(
     return {read: int(label) for read, label in zip(read_names, labels, strict=True)}
 
 
-def pairwise_alignment_lengths(ava_alignments:list[pysam.AlignedSegment], read_names:list[str]) -> np.ndarray:
+def pairwise_alignment_lengths(
+    ava_alignments: list[pysam.AlignedSegment], read_names: list[str]
+) -> np.ndarray:
     """Calculate pairwise alignment lengths from all-vs-all alignments.
 
     Args:
@@ -1472,7 +1485,7 @@ def consensus_while_clustering(
                 _pairwise_aln_lengths = pairwise_alignment_lengths(
                     ava_alignments=all_vs_all_alignments, read_names=sim_read_names
                 )
-                
+
                 _representative = (
                     find_representative_read(
                         similarity_matrix=similarity_matrix,
@@ -3379,4 +3392,3 @@ def main():
 if __name__ == "__main__":
     main()
 # %%
-

src/svirlpool/signalprocessing/signaldepths_to_signalstrength.py

@@ -67,7 +67,11 @@ def genomic_distance_score(d: int, flatness: float = 15.0) -> float:
 
 
 def delta(
-    subj: npt.NDArray[np.int32], obj: npt.NDArray[np.int32], flatness: float
+    subj: npt.NDArray[np.int32],
+    obj: npt.NDArray[np.int32],
+    flatness: float,
+    bnd_proximity: int = 100,
+    bnd_repeat_proximity: int = 500,
 ) -> float:
     """returns 1 if no distance between two loci, else x >= 0"""
     distance = np.uint32(2**31 - 1)
@@ -85,9 +89,9 @@ def delta(
         case (0, 2):  # INS vs DELR
             pass
         case (0, 3):  # INS vs BNDL
-            distance = distance_collapsed_sv_signals(subj, obj)
+            pass
         case (0, 4):  # INS vs BNDR
-            distance = distance_collapsed_sv_signals(subj, obj)
+            pass
         case (1, 0):  # DELL vs INS
             pass
         case (1, 1):  # DELL vs DELL
@@ -98,7 +102,7 @@ def delta(
         case (1, 3):  # DELL vs BNDL
             pass
         case (1, 4):  # DELL vs BNDR
-            distance = distance_collapsed_sv_signals(subj, obj)
+            pass
         case (2, 0):  # DELR vs INS
             pass
         case (2, 1):  # DELR vs DELL
@@ -107,7 +111,7 @@ def delta(
             distance = distance_collapsed_sv_signals(subj, obj)
             size_score = size_similarity_score(subj, obj)
         case (2, 3):  # DELR vs BNDL
-            distance = distance_collapsed_sv_signals(subj, obj)
+            pass
         case (2, 4):  # DELR vs BNDR
             pass
         case (3, 0):  # BNDL vs INS
@@ -118,6 +122,9 @@ def delta(
             pass
         case (3, 3):  # BNDL vs BNDL
             distance = distance_collapsed_sv_signals(subj, obj)
+            _in_same_repeat = subj[-2] == obj[-2] and int(subj[-2]) >= 0
+            if distance > (bnd_repeat_proximity if _in_same_repeat else bnd_proximity):
+                return 0.0
         case (3, 4):  # BNDL vs BNDR
             pass
         case (4, 0):  # BNDR vs INS
@@ -130,12 +137,15 @@ def delta(
             pass
         case (4, 4):  # BNDR vs BNDR
             distance = distance_collapsed_sv_signals(subj, obj)
+            _in_same_repeat = subj[-2] == obj[-2] and int(subj[-2]) >= 0
+            if distance > (bnd_repeat_proximity if _in_same_repeat else bnd_proximity):
+                return 0.0
     if distance < 0:
         return 1.0
     elif distance == 0:
         return 1.0
     else:
-        distance_score = genomic_distance_score(distance, flatness)
+        distance_score = genomic_distance_score(int(distance), flatness)
         return max(distance_score, size_score)
 
 
@@ -146,6 +156,8 @@ def pseudo_convolve(
     warning_collapsed_region_size: int = 2_000,
     high_density_threshold: int = 200,
     high_density_score: float = 10.0,
+    bnd_proximity: int = 100,
+    bnd_repeat_proximity: int = 500,
 ) -> np.ndarray:
     """convolves over a numpy array of signals and returns the sum of weighted distances
     by a kernel function (radius) which is dynamically expanded if it overlaps any repeats.
@@ -207,7 +219,13 @@ def pseudo_convolve(
         for l in range(j, k):  # noqa: E741
             if l == i:
                 continue
-            next_val = delta(subj=np_signals[i], obj=np_signals[l], flatness=flatness)
+            next_val = delta(
+                subj=np_signals[i],
+                obj=np_signals[l],
+                flatness=flatness,
+                bnd_proximity=bnd_proximity,
+                bnd_repeat_proximity=bnd_repeat_proximity,
+            )
             values[i] += next_val
 
         i += 1
@@ -413,6 +431,8 @@ def signaldepths_to_signalstrength(
     kernel_signal_radius: int,
     slop_repeats_radius: int,
     tmp_dir_path: Path | None = None,
+    bnd_proximity: int = 100,
+    bnd_repeat_proximity: int = 500,
 ) -> None:
     csv.field_size_limit(sys.maxsize)
     # --- load and prepare data --- #
@@ -520,6 +540,8 @@ def signaldepths_to_signalstrength(
                 np_signals=np_signals_chr,
                 kernelsize=kernel_signal_radius,
                 flatness=flatness_distance,
+                bnd_proximity=bnd_proximity,
+                bnd_repeat_proximity=bnd_repeat_proximity,
             )
 
             # Write results immediately
@@ -573,6 +595,8 @@ def run(args, **kwargs):
         kernel_signal_radius=args.kernel_signal_radius,
         threads=args.threads,
         slop_repeats_radius=args.slop_repeats_radius,
+        bnd_proximity=args.bnd_proximity,
+        bnd_repeat_proximity=args.bnd_repeat_proximity,
     )
 
 
@@ -635,6 +659,20 @@ def get_parser():
         default=50,
         help="Radius of slop for repeats.",
     )
+    parser.add_argument(
+        "--bnd-proximity",
+        type=int,
+        required=False,
+        default=100,
+        help="Maximum distance (bp) between two BND signals for them to score against each other.",
+    )
+    parser.add_argument(
+        "--bnd-repeat-proximity",
+        type=int,
+        required=False,
+        default=500,
+        help="Maximum distance (bp) between two BND signals within the same tandem repeat for them to score against each other.",
+    )
     return parser