Severus 1.5 gives IndexError: list index out of range

[byoo]

Hello!

Could you help with the error below? I ran Severus 1.5 on a single sample/tumor sequenced by PacBio as follows:
/sched/software/python/conda/envs/severus_150/bin/severus --target-bam caVal-093.GRCh38.haplotagged.bam --out-dir caVal-093 -t 16 --min-support 1 --phasing-vcf caVal-093.GRCh38.deepvariant.phased.vcf.gz --vntr-bed human_GRCh38_no_alt_analysis_set.trf.bed

Below is the log:

[2025-06-26 22:47:28] INFO: Starting Severus 1.5
[2025-06-26 22:47:55] INFO: Parsing reads from caVal-093.GRCh38.haplotagged.bam
[2025-06-26 22:50:40] INFO: Total read length: 77398671581
[2025-06-26 22:50:40] INFO: Total aligned length: 74747628317 (0.97)
[2025-06-26 22:50:40] INFO: Read N50 / N90: 11112 / 5307
[2025-06-26 22:50:40] INFO: Alignments N50 / N90: 11047 / 5265
[2025-06-26 22:50:40] INFO: Read error rate (Q25 / Q50 / Q75): 0.0000 / 0.0020 / 0.0040
[2025-06-26 22:50:40] INFO: Read mismatch rate (Q25 / Q50 / Q75): 0.0000 / 0.0000 / 0.0010
[2025-06-26 22:51:27] INFO: Computing read quality
[2025-06-26 22:51:40] INFO: Annotating reads
[2025-06-26 22:51:48] INFO: Computing coverage histogram
[2025-06-26 22:51:49] INFO: Median coverage by PASS reads for caVal-093.GRCh38.haplotagged.bam (H1 / H2 / H0): 10.0 / 10.0 / 0.0
[2025-06-26 22:51:49] INFO: Extracting split alignments
[2025-06-26 22:51:52] INFO: Extracting clipped reads
[2025-06-26 22:51:52] INFO: Starting breakpoint detection
[2025-06-26 22:52:06] INFO: Clustering unmapped insertions
/sched/software/python/conda/envs/severus_150/lib/python3.12/site-packages/numpy/_core/fromnumeric.py:3860: RuntimeWarning: Mean of empty slice.
return _methods._mean(a, axis=axis, dtype=dtype,
/sched/software/python/conda/envs/severus_150/lib/python3.12/site-packages/numpy/_core/_methods.py:145: RuntimeWarning: invalid value encountered in scalar divide
ret = ret.dtype.type(ret / rcount)
[2025-06-26 22:52:29] INFO: Starting compute_bp_coverage
[2025-06-26 22:52:49] INFO: Filtering breakpoints
[2025-06-26 22:53:00] INFO: Writing breakpoints
[2025-06-26 22:53:06] INFO: Preparing outputs for all_SVs
[2025-06-26 22:53:06] INFO: Computing segment coverage
[2025-06-26 22:54:30] INFO: Total phased length: 2176060957
[2025-06-26 22:54:30] INFO: Phase blocks N50: 330693
Traceback (most recent call last):
File "/sched/software/python/conda/envs/severus_150/bin/severus", line 33, in
sys.exit(load_entry_point('severus==1.5', 'console_scripts', 'severus')())
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/sched/software/python/conda/envs/severus_150/lib/python3.12/site-packages/severus/main.py", line 260, in main
output_graphs(double_breaks, coverage_histograms, thread_pool, target_genomes, control_genomes, genome_ids, ref_lengths, args)
File "/sched/software/python/conda/envs/severus_150/lib/python3.12/site-packages/severus/build_graph.py", line 655, in output_graphs
(genomic_segments, adj_segments) = get_genomic_segments(double_breaks, coverage_histograms, args.phase_vcf,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/sched/software/python/conda/envs/severus_150/lib/python3.12/site-packages/severus/breakpoint_finder.py", line 1741, in get_genomic_segments
calc_gen_segments(d_breaks, coverage_histograms,ref_lengths, min_ref_flank, max_genomic_length, db_segments, ref_adj)
File "/sched/software/python/conda/envs/severus_150/lib/python3.12/site-packages/severus/breakpoint_finder.py", line 1858, in calc_gen_segments
if db.bp_2.position == pos_ls[ind]:
~~~~~~^^^^^
IndexError: list index out of range

[aysegokce]

Hello Byunggil,

Could you try the current version in ayse_dev?

Best
Ayse

[byoo]

Hi Ayse,

It failed with the ayse_dev branch as well. Here's the log:

[2025-07-10 17:02:43] INFO: Starting Severus 1.6
[2025-07-10 17:04:09] INFO: Parsing reads from caVal-093.GRCh38.haplotagged.bam
[2025-07-10 17:08:23] INFO:     Total read length: 77398671581
[2025-07-10 17:08:23] INFO:     Total aligned length: 74747628317 (0.97)
[2025-07-10 17:08:23] INFO:     Read N50 / N90: 11112 / 5307
[2025-07-10 17:08:23] INFO:     Alignments N50 / N90: 11047 / 5265
[2025-07-10 17:08:23] INFO:     Read error rate (Q25 / Q50 / Q75): 0.0000 / 0.0020 / 0.0040
[2025-07-10 17:08:23] INFO:     Read mismatch rate (Q25 / Q50 / Q75): 0.0000 / 0.0000 / 0.0010
[2025-07-10 17:09:01] INFO: Computing read quality
[2025-07-10 17:09:12] INFO: Annotating reads
[2025-07-10 17:09:19] INFO: Computing coverage histogram
[2025-07-10 17:09:20] INFO:     Median coverage by PASS reads for caVal-093.GRCh38.haplotagged.bam (H1 / H2 / H0): 10.0 / 10.0 / 0.0
[2025-07-10 17:09:20] INFO: Extracting split alignments
[2025-07-10 17:09:20] INFO: Extracting clipped reads
[2025-07-10 17:09:23] INFO: Starting breakpoint detection
[2025-07-10 17:09:34] INFO: Clustering unmapped insertions
/sched/software/python/conda/envs/severus_dev/lib/python3.12/site-packages/numpy/_core/fromnumeric.py:3860: RuntimeWarning: Mean of empty slice.
  return _methods._mean(a, axis=axis, dtype=dtype,
/sched/software/python/conda/envs/severus_dev/lib/python3.12/site-packages/numpy/_core/_methods.py:144: RuntimeWarning: invalid value encountered in scalar divide
  ret = ret.dtype.type(ret / rcount)
[2025-07-10 17:09:53] INFO: Starting compute_bp_coverage
[2025-07-10 17:10:16] INFO: Filtering breakpoints
[2025-07-10 17:10:25] INFO: Writing breakpoints
[2025-07-10 17:10:31] INFO: Preparing outputs for all_SVs
[2025-07-10 17:10:31] INFO:     Computing segment coverage
[2025-07-10 17:12:24] INFO:     Total phased length: 2176060957
[2025-07-10 17:12:24] INFO:     Phase blocks N50: 330693
Traceback (most recent call last):
  File "/sched/software/python/conda/envs/severus_dev/bin/severus", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/sched/software/python/conda/envs/severus_dev/lib/python3.12/site-packages/severus/main.py", line 263, in main
    output_graphs(double_breaks, coverage_histograms, thread_pool, target_genomes, control_genomes, genome_ids, ref_lengths, args)
  File "/sched/software/python/conda/envs/severus_dev/lib/python3.12/site-packages/severus/build_graph.py", line 655, in output_graphs
    (genomic_segments, adj_segments) = get_genomic_segments(double_breaks, coverage_histograms, args.phase_vcf,
                                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/sched/software/python/conda/envs/severus_dev/lib/python3.12/site-packages/severus/breakpoint_finder.py", line 1786, in get_genomic_segments
    calc_gen_segments(d_breaks, coverage_histograms,ref_lengths, max_genomic_length, db_segments, ref_adj)
  File "/sched/software/python/conda/envs/severus_dev/lib/python3.12/site-packages/severus/breakpoint_finder.py", line 1903, in calc_gen_segments
    if db.bp_2.position == pos_ls[ind]:
                           ~~~~~~^^^^^
IndexError: list index out of range

[byoo]

FYI, the data can be shared if that makes the troubleshooting easier. If you help me locate the region causing the problem, I will subset the BAM/VCF and send you quickly.

[aysegokce]

That would be really helpful! Thank you!
Ayse

[byoo]

Thank you! The whole data would take a few days. Is it possible to find the region that cause the problem? This will allow me to send the data via email quickly.

I tried to print to find out the regions of interest to subset but I wasn't sure how to do that.
But I see that ind and len(pos_ls) being equal (1897) is problematic; it appears that Severus assumes that db.bp_2.position is in pos_ls but it's not the case here and led to the index out of range error.

If you like to have the whole BAM/VCF, I will initiate the transfer.

ind: 1897, len(pos_ls): 1897, db_cov: 1, bp2_len: 3313.0
ref_name: chr13 db.bp_1.position: 113669568 db.bp_2.position: 145047879
Traceback (most recent call last):
  File "/sched/software/python/conda/envs/severus_dev/bin/severus", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/sched/software/python/conda/envs/severus_dev/lib/python3.12/site-packages/severus/main.py", line 263, in main
    output_graphs(double_breaks, coverage_histograms, thread_pool, target_genomes, control_genomes, genome_ids, ref_lengths, args)
  File "/sched/software/python/conda/envs/severus_dev/lib/python3.12/site-packages/severus/build_graph.py", line 655, in output_graphs
    (genomic_segments, adj_segments) = get_genomic_segments(double_breaks, coverage_histograms, args.phase_vcf,
                                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/sched/software/python/conda/envs/severus_dev/lib/python3.12/site-packages/severus/breakpoint_finder.py", line 1786, in get_genomic_segments
    calc_gen_segments(d_breaks, coverage_histograms,ref_lengths, max_genomic_length, db_segments, ref_adj)
  File "/sched/software/python/conda/envs/severus_dev/lib/python3.12/site-packages/severus/breakpoint_finder.py", line 1905, in calc_gen_segments
    if db.bp_2.position == pos_ls[ind]:
                           ~~~~~~^^^^^
IndexError: list index out of range

[aysegokce]

Could you add following code before the 'if db.bp_2.position == pos_ls[ind]:'?

if ind == len(pos_ls):
     print((ref_name, db.bp_2.ref_id, db.bp_2.position,db.bp_1.ref_id, db.bp_1.position, pos_ls[-1]))

Best
Ayse

[byoo]

It prints the following info:
('chr13', 'chr13', 145047879, 'chr8', 113669568, 114365328)

[aysegokce]

Thank you!
Is it possible to share a BAM with chr8:135138636-145138636 and chr13:104364328-114364328 and the breakpoints_double.csv from you runs?

Thank you!