Changeset 5775

Timestamp:

Dec 6, 2019, 7:56:14 AM (3 years ago)

Author:

Nicklas Nordborg

Message:

References #1210: Variant calling fails if no variants are found

Added a test to the script after VarDict that checks the number of variants. If there are no variants in the VCF, the rest of the pipeline is skipped and the raw file is simply copied to the expected places.

File:

• extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/grid/VariantCallingJobCreator.java (modified) (7 diffs)

extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/grid/VariantCallingJobCreator.java

@@ -357 +357 @@
         script.comment("Skipping raw variant calling since we can use existing file");
         script.cmd("zcat ${BamFolder}/" + vcfName + " > tmp/variants-raw-2.vcf");
+        // 'grep' to find number of variants The ' || true' is critical since 'grep' will exit with code 1 if there are no matches
+        script.cmd("NUMVARIANTS=`grep -c -v '^#' tmp/variants-raw-2.vcf || true`");
         script.newLine();
       }
@@ -395 +397 @@
         script.progress(30, "Running VarDict (${NumThreads} threads)");
         script.cmd(vardictCmd);
-        script.newLine();
-      
+        // 'grep' to find number of variants The ' || true' is critical since 'grep' will exit with code 1 if there are no matches
+        script.cmd("NUMVARIANTS=`grep -c -v '^#' tmp/variants-raw-1.vcf || true`");
+        script.newLine();
+      
+        script.cmd("if [ $NUMVARIANTS -gt 0 ]; then");
         script.comment("Calculate GC content and DNA strength using custom script");
         script.comment("We need to get the reference sequence within 50 base-pairs of each entry in the VCF");
@@ -406 +411 @@
         script.cmd("tabix tmp/gc-50.vcf.gz");
         script.cmd("./stderrwrap.sh " + vcfanno_path + " gc_stat.toml tmp/variants-raw-1.vcf > tmp/variants-raw-2.vcf 3>> vcfanno.out");
+        script.cmd("else");
+        script.progress(33, "Before copy");
+        script.comment("No variants found, simply copy the raw file");
+        script.cmd("cp tmp/variants-raw-1.vcf tmp/variants-raw-2.vcf");
+        script.progress(34, "After copy");
+        script.cmd("fi");
         script.cmd("cat tmp/variants-raw-2.vcf | bgzip -c > resultsraw/variants-raw.vcf.gz");
         script.newLine();
@@ -412 +423 @@
       if (!rawOnly)
       {
+        script.cmd("if [ $NUMVARIANTS -gt 0 ]; then");
         script.comment("Calculate distance to closest G5 SNP and add back");
         script.comment("annotations DIST_SNP and DIST_DIV to the variants VCF");
@@ -456 +468 @@
         script.cmd(snpSiftCmd);
         script.newLine();
+        script.cmd("else");
+        script.cmd("  cat tmp/variants-raw-2.vcf | bgzip -c > resultsfilter/variants-annotated.vcf.gz");
+        script.cmd("  cp tmp/variants-raw-2.vcf resultsfilter/variants-filtered.vcf");
+        script.cmd("fi");
       }
 
@@ -463 +479 @@
         script.cmd("\\cp resultsraw/variants-* ${BamFolder}");  // "\cp" will enable overwriting existing files
         script.cmd("ls -1 resultsraw/variants-* >> ${WD}/files.out");
-        script.cmd("echo \"Callable bases: `awk -F'\t' 'BEGIN{SUM=0}{ SUM+=$3-$2 } END{print SUM}' resultsraw/variants-callable.bed`\" >> ${WD}/stats.out");
-        script.cmd("echo \"Raw variants: `zcat resultsraw/variants-raw.vcf.gz | grep -v '^#\' | wc -l`\" >> ${WD}/stats.out");
+        script.cmd("echo \"Callable bases: `awk -F'\\t' 'BEGIN{SUM=0}{ SUM+=$3-$2 } END{print SUM}' resultsraw/variants-callable.bed`\" >> ${WD}/stats.out");
+        script.cmd("echo \"Raw variants: `zcat resultsraw/variants-raw.vcf.gz | grep -c -v '^#'`\" >> ${WD}/stats.out");
         if (externalGroup != null)
         {
@@ -476 +492 @@
         script.cmd("cp resultsfilter/* ${FilteredFolder}");
         script.cmd("ls -1 ${FilteredFolder}/* >> ${WD}/files.out");
-        script.cmd("echo \"Annotated variants: `zcat resultsfilter/variants-annotated.vcf.gz | grep -v '^#\' | wc -l`\" >> ${WD}/stats.out");
-        script.cmd("echo \"Filtered variants: `cat resultsfilter/variants-filtered.vcf | grep -v '^#\' | wc -l`\" >> ${WD}/stats.out");
+        script.cmd("echo \"Annotated variants: `zcat resultsfilter/variants-annotated.vcf.gz | grep -c -v '^#'`\" >> ${WD}/stats.out");
+        script.cmd("echo \"Filtered variants: `cat resultsfilter/variants-filtered.vcf | grep -c -v '^#'`\" >> ${WD}/stats.out");
         if (externalGroup != null)
         {