Updating filter vcf script to handle edge cases

mgymrek · Jul 1, 2016 · 47da4ca · 47da4ca
1 parent eead65a
commit 47da4ca
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 6 deletions.
diff --git 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
@@ -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"