From c84c66cc5bb388262ee1eda181feac94449eef3f Mon Sep 17 00:00:00 2001
From: Venkat Malladi <venkat.malladi@utsouthwestern.edu>
Date: Wed, 12 Feb 2020 00:25:08 -0600
Subject: [PATCH] Solve the minimum fragment length to be at least 50bp or
 through an error.

---
 .gitignore                             |  1 -
 test_data/test_cross.qc                |  1 +
 workflow/scripts/call_peaks_macs.py    | 18 ++++++++++++++----
 workflow/tests/test_call_peaks_macs.py | 19 +++++++++++++++++++
 4 files changed, 34 insertions(+), 5 deletions(-)
 create mode 100644 test_data/test_cross.qc

diff --git a/.gitignore b/.gitignore
index bf385c6..7756350 100644
--- a/.gitignore
+++ b/.gitignore
@@ -108,5 +108,4 @@ report*.html*
 timeline*.html*
 /workflow/output/*
 /work/*
-/test_data/*
 /.nextflow/*
diff --git a/test_data/test_cross.qc b/test_data/test_cross.qc
new file mode 100644
index 0000000..6fbc0ad
--- /dev/null
+++ b/test_data/test_cross.qc
@@ -0,0 +1 @@
+Test.20.tagAlign.gz	18588987	0,20,33	0.211525291335199,0.211232019956852,0.211139666755398	35	0.2123067	1500	0.209429	1.01001	0.7284536       0
diff --git a/workflow/scripts/call_peaks_macs.py b/workflow/scripts/call_peaks_macs.py
index 6981262..f7527eb 100644
--- a/workflow/scripts/call_peaks_macs.py
+++ b/workflow/scripts/call_peaks_macs.py
@@ -138,10 +138,20 @@ def call_peaks_macs(experiment, xcor, control, prefix, genome_size, chrom_sizes)
     with open(xcor, 'r') as xcor_fh:
         firstline = xcor_fh.readline()
         frag_lengths = firstline.split()[2]  # third column
-        fragment_length = frag_lengths.split(',')[0]  # grab first value
-        logger.info("Fraglen %s", fragment_length)
-        if int(fragment_length) < 1:
-            fragment_length = "200"
+        frag_lengths_array = frag_lengths.split(',')
+        fragment_length = 0
+        fragment = False
+        # Loop through all values of fragment length
+        for f in frag_lengths.split(','):
+            fragment_length = f
+            logger.info("Fraglen %s", fragment_length)
+            if int(fragment_length) > 50:
+                fragment = True
+                break
+
+        if fragment == False:
+            logger.info('Error in cross-correlation analysis: %s', frag_lengths_array)
+            raise Exception("Error in cross-correlation analysis: %s" % frag_lengths_array)
 
     # Generate narrow peaks and preliminary signal tracks
 
diff --git a/workflow/tests/test_call_peaks_macs.py b/workflow/tests/test_call_peaks_macs.py
index cd94e17..3b4404d 100644
--- a/workflow/tests/test_call_peaks_macs.py
+++ b/workflow/tests/test_call_peaks_macs.py
@@ -3,10 +3,29 @@
 import pytest
 import os
 import utils
+from io import StringIO
+import call_peaks_macs
 
 test_output_path = os.path.dirname(os.path.abspath(__file__)) + \
                 '/../output/callPeaksMACS/'
 
+test_data_path = os.path.dirname(os.path.abspath(__file__)) + \
+		'/../../test_data/'
+        
+
+@pytest.mark.unit
+def test_fragment_length():
+    experiment = "experiment.tagAlign.gz"
+    control = "control.tagAlign.gz"
+    prefix = 'test'
+    genome_size = 'hs'
+    chrom_sizes = 'genomefile.txt'
+    cross_qc = os.path.join(test_data_path, 'test_cross.qc')
+    with pytest.raises(Exception) as excinfo:
+        call_peaks_macs.call_peaks_macs(experiment, cross_qc, control, prefix, genome_size, chrom_sizes)
+    assert str(excinfo.value) == "Error in cross-correlation analysis: ['0', '20', '33']"
+
+
 
 @pytest.mark.singleend
 def test_fc_signal_singleend():
-- 
GitLab