From 47da4cac554952f14a78bb4230d39efa1eaf8b0b Mon Sep 17 00:00:00 2001
From: mgymrek <mgymrek@mit.edu>
Date: Fri, 1 Jul 2016 14:57:21 -0400
Subject: [PATCH] Updating filter vcf script to handle edge cases

---
 scripts/lobSTR_filter_vcf.py | 16 ++++++++++++----
 src/main_allelotype.cpp      |  3 +--
 2 files changed, 13 insertions(+), 6 deletions(-)

diff --git a/scripts/lobSTR_filter_vcf.py b/scripts/lobSTR_filter_vcf.py
index 4537206..d5897d5 100755
--- a/scripts/lobSTR_filter_vcf.py
+++ b/scripts/lobSTR_filter_vcf.py
@@ -128,6 +128,11 @@ def GetWriter(reader, filters):
     numallcalls = 0
     numpasscalls = 0
 
+    # Test if there are any samples left
+    if len(reader.samples)-len(ignore_samples) == 0:
+        MSG("Exiting, no samples remaining")
+        sys.exit(0)
+
     # Go through each record
     for record in reader:
         # Set up new sample info. only need to do this once per record
@@ -185,7 +190,6 @@ def GetWriter(reader, filters):
             str_length = record.INFO["END"]-record.POS+1
             if str_length > filters[KEY_LOCMAXREFLEN]["Value"]:
                 record.add_filter(KEY_LOCMAXREFLEN+str(filters[KEY_LOCMAXREFLEN]["Value"]))
-                continue
         if KEY_LOCCALLRATE in filters:
             if len(ignore_samples) > 0:
                 call_rate = len(coverages)*1.0/(len(reader.samples)-len(ignore_samples))
@@ -205,6 +209,13 @@ def GetWriter(reader, filters):
         allfilters.append(":".join(record.FILTER))
         writer.write_record(record)
         
+    if len(allfilters) == 0:
+        MSG("No loci detected")
+        sys.exit(0)
+    if numallcalls == 0:
+        MSG("No calls detected")
+        sys.exit(0)
+        
     MSG("### Locus Filter counts ###")
     for f in set(allfilters):
         MSG("%s: %s (%s)"%(f, allfilters.count(f), allfilters.count(f)*1.0/len(allfilters)))
@@ -213,6 +224,3 @@ def GetWriter(reader, filters):
     MSG("  Mean coverage: %s"%np.mean(allcoverages))
     MSG("Number of calls passing filtering: %s (%s)"%(numpasscalls, numpasscalls*1.0/numallcalls))
     MSG("  Mean coverage: %s"%np.mean(passcoverages))
-        
-        
-
diff --git a/src/main_allelotype.cpp b/src/main_allelotype.cpp
index ecb9da9..2a5e9d8 100644
--- a/src/main_allelotype.cpp
+++ b/src/main_allelotype.cpp
@@ -90,8 +90,7 @@ void show_help() {
 	   << "--no-rmdup:                     don't remove pcr duplicates before allelotyping.\n"
      << "--realign:                      Redo local realignment. Useful if using alignments generated\n"
      << "                                by other tools.\n"
-	   << "--min-het-freq <FLOAT>:         minimum frequency to make a heterozygous call\n"
-	   << "                                (default: 0.1)\n"
+     << "--min-het-freq <FLOAT>:         minimum frequency to make a heterozygous call. (Default: " << min_het_freq << ")\n"
 	   << "--haploid <chrX,[chrY,...]>:    comma-separated list of chromosomes\n"
 	   << "                                that should be forced to have homozygous\n"
 	   << "                                calls. Specify --haploid all if the organism\n"