Fix error when too few reads on chromosome

src/naibr/distributions.py

@@ -1,5 +1,6 @@
 import collections
 import itertools
+import logging
 import os
 
 import matplotlib as mpl
@@ -15,6 +16,8 @@
 
 mpl.use("Agg")
 
+logger = logging.getLogger(__name__)
+
 
 class NegBin:
     def __init__(self, p=0.1, r=10):
@@ -152,6 +155,7 @@ def get_linkedread_distributions(linkedreads_by_barcode, configs):
             get_pairwise_overlap(barcode_linkedreads, barcode_overlap, configs)
 
     if len(linkedreads) < 100:
+        logger.warning("Too few linked reads to estimate distributions")
         return None, None, None
 
     p_rate = get_rate_distr(linkedreads)

src/naibr/get_reads.py

@@ -248,7 +248,7 @@ def get_candidates_from_discs(discs, barcode_overlap, configs):
 def get_candidates(discs, reads_by_barcode, configs):
     p_len, p_rate, barcode_overlap, barcode_linkedreads = get_distributions(reads_by_barcode, configs)
     if p_len is None or p_rate is None:
-        return [], None, None, None
+        return [], None, None, collections.defaultdict(list)
 
     candidates = get_candidates_from_discs(discs, barcode_overlap, configs)
     return candidates, p_len, p_rate, barcode_linkedreads

tests/data/make_testdata.sh

@@ -41,3 +41,8 @@ bam="https://s3-us-west-2.amazonaws.com/10x.files/samples/genome/HCC1954T_WGS_21
 region="chr3:101,508,502-101,725,227"
 samtools view "${bam}" "${region}" -o HCC1954T_10xGenomics_chr3_DUP.bam
 samtools index HCC1954T_10xGenomics_chr3_DUP.bam
+
+
+# Subsample BAM to contain only a single barcode
+samtools view -d BX:AACGGTTCACGAAACG-1 tests/data/HCC1954T_10xGenomics_chr21_INV.bam -bh -o tests/data/HCC1954T_10xGenomics_chr21_INV.one_barcode.bam
+samtools index tests/data/HCC1954T_10xGenomics_chr21_INV.one_barcode.bam

tests/test_run.py

@@ -13,6 +13,7 @@
 BAM_DEL = Path("tests/data/NA24385_10xGenomics_chr1_DEL.bam").absolute()
 BAM_DUP = Path("tests/data/HCC1954T_10xGenomics_chr3_DUP.bam").absolute()
 BAM_TRA = Path("tests/data/HCC1954T_10xGenomics_chr_12_20_TRA.bam").absolute()
+BAM_SMALL = Path("tests/data/HCC1954T_10xGenomics_chr21_INV.one_barcode.bam").absolute()
 
 
 def same_pairwise_elements(list1, list2):
@@ -231,3 +232,15 @@ def test_call_interchromosomal(tmp_path):
         assert nas_pass[0].chrm2 == "chr20"
         assert abs(nas_pass[0].break1 - 40393752) < 500
         assert abs(nas_pass[0].break2 - 53974890) < 500
+
+
+def test_too_few_barcodes(tmp_path):
+    config_file = io.StringIO(f"bam_file={BAM_SMALL}\noutdir={tmp_path}\nDEBUG=True\n")
+    configs = Configs.from_file(config_file)
+
+    exitcode = run(configs)
+    assert exitcode == 0
+    output = tmp_path / "NAIBR_SVs.bedpe"
+    with open(output) as bedpe_reader:
+        nas = list(iter_novel_adjacencies(bedpe_reader))
+        assert len(nas) == 0