Optionally emit SV lengths

setup.py

@@ -2,7 +2,7 @@
 
 setup(
  name='sniffles',
- version='2.4',
+ version='2.4.1',
  packages=find_packages(),
  url='https:/fritzsedlazeck/Sniffles',
  license='MIT',
@@ -11,4 +11,4 @@
  description='A fast structural variation caller for long-read sequencing data',
  long_description=open('README.md').read(),
  long_description_content_type='text/markdown',
-)
+)

src/sniffles/config.py

@@ -3,7 +3,7 @@
 # Sniffles2
 # A fast structural variant caller for long-read sequencing data
 #
-# Created:  18.10.2021
+# Created: 18.10.2021
 # Author: Moritz Smolka
 # Maintainer: Hermann Romanek
 # Contact: [email protected]
@@ -13,15 +13,14 @@
 import sys
 import datetime
 import argparse
-from collections import defaultdict
 
 from typing import Union, Optional
 
 from sniffles import util
 from sniffles.region import Region
 
 VERSION = "Sniffles2"
-BUILD = "2.4"
+BUILD = "2.4.1"
 SNF_VERSION = "S2_rc4"
 
 
@@ -224,6 +223,8 @@ def add_mosaic_args(self, parser):
 
  def add_developer_args(self, parser):
  developer_args = parser.add_argument_group("Developer parameters")
+
+ developer_args.add_argument("--dev-emit-sv-lengths", default=False, action="store_true", help=argparse.SUPPRESS)
  developer_args.add_argument("--dev-cache", default=False, action="store_true", help=argparse.SUPPRESS)
  developer_args.add_argument("--dev-cache-dir", metavar="PATH", type=str, default=None, help=argparse.SUPPRESS)
  developer_args.add_argument("--dev-debug-svtyping", default=False, action="store_true", help=argparse.SUPPRESS)

src/sniffles/sv.py

@@ -49,6 +49,7 @@ class SVCall:
 
  svtype: str
  svlen: int
+ svlens: list[int]
  end: int
  genotypes: dict[int, tuple]
 
@@ -268,6 +269,9 @@ def call(self, config, task) -> Optional[SVCall]:
 
  svcall_pos = int(util.median(cand.pos for cand in self.candidates))
  svcall_svlen = int(util.median(cand.svlen for cand in self.candidates))
+ svcall_svlens: list[int] = [length in cand_lengths
+ for cand_lengths in self.candidates
+ for length in cand_lengths.svlens]
  svcall_alt = first_cand.alt
  svcall_alt_mindist = abs(len(svcall_alt) - svcall_svlen)
  if first_cand.svtype == "INS":
@@ -290,6 +294,7 @@ def call(self, config, task) -> Optional[SVCall]:
  info=dict(),
  svtype=first_cand.svtype,
  svlen=svcall_svlen if config.dev_combine_medians else first_cand.svlen,
+ svlens=svcall_svlens,
  end=svcall_end if config.dev_combine_medians else first_cand.end,
  genotypes=genotypes,
  precise=sum(int(cand.precise) for cand in self.candidates) / float(len(self.candidates)) > 0.5,
@@ -338,7 +343,7 @@ def call_from(cluster, config, keep_qc_fails, task):
  qc = True
 
  svlen = util.center(v.svlen for v in leads)
-
+ svlens = [v.svlen for v in leads]
  if abs(svlen) < config.minsvlen_screen:
  return
 
@@ -394,6 +399,7 @@ def call_from(cluster, config, keep_qc_fails, task):
  info=dict(),
  svtype=svtype,
  svlen=svlen,
+ svlens=svlens,
  end=svend,
  genotypes=dict(),
  precise=precise,

src/sniffles/vcf.py

@@ -81,6 +81,9 @@ def __init__(self, config: SnifflesConfig, handle):
  if config.qc_nm_measure:
  self.info_order.append("NM")
 
+ if config.dev_emit_sv_lengths:
+ self.info_order.append("SVLENGTHS")
+
  self.default_genotype = config.genotype_none
 
  # Add phasing if needed
@@ -151,6 +154,7 @@ def write_header(self, contigs_lengths):
  self.write_header_line('INFO=<ID=IMPRECISE,Number=0,Type=Flag,Description="Structural variation with imprecise breakpoints">')
  self.write_header_line('INFO=<ID=MOSAIC,Number=0,Type=Flag,Description="Structural variation classified as putative mosaic">')
  self.write_header_line('INFO=<ID=SVLEN,Number=1,Type=Integer,Description="Length of structural variation">')
+ self.write_header_line('INFO=<ID=SVLENGTH,Number=.,Type=Integer,Description="Lengths of structural variation (all)">')
  self.write_header_line('INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variation">')
  self.write_header_line('INFO=<ID=CHR2,Number=1,Type=String,Description="Mate chromsome for BND SVs">')
  self.write_header_line('INFO=<ID=SUPPORT,Number=1,Type=Integer,Description="Number of reads supporting the structural variation">')
@@ -218,6 +222,7 @@ def write_call(self, call):
  infos = {
  "SVTYPE": call.svtype,
  "SVLEN": call.svlen,
+ "SVLENGTHS": ",".join(map(str, call.svlens)),
  "END": end,
  "SUPPORT": call.support,
  "RNAMES": call.rnames if self.config.output_rnames else None,
@@ -229,6 +234,7 @@ def write_call(self, call):
 
  if call.svtype == "BND":
  infos["SVLEN"] = None
+ infos["SVLENGTHS"] = None
  infos["END"] = None
 
  infos_ordered = ["PRECISE" if call.precise else "IMPRECISE"]
@@ -343,6 +349,9 @@ def read_svs_iter(self):
 
  if "SVLEN" in info_dict:
  call.svlen = int(info_dict["SVLEN"])
+ if "SVLENGTHS`" in info_dict:
+ call.svlens = info_dict["SVLENGTHS"]
+
  if "END" in info_dict:
  call.end = int(info_dict["END"])