report pct_top_hit_of_total_reads as one of the outputs of align_and_…

…count report pct_top_hit_of_total_reads as one of the outputs of align_and_count, calculated as: [(mapped read count for the spike-in with the highest count) / (total number of reads in the bam file)]
broadinstitute · Jan 23, 2025 · df03eef · df03eef
1 parent 397f131
commit df03eef
Show file tree

Hide file tree

Showing 2 changed files with 14 additions and 4 deletions.
diff --git a/pipes/WDL/tasks/tasks_reports.wdl b/pipes/WDL/tasks/tasks_reports.wdl
@@ -480,6 +480,14 @@ task align_and_count {
     else
       PCT_OF_INPUT_READS_MAPPED=$( echo "scale=3; 100 * ($TOTAL_COUNT_OF_LESSER_HITS + $TOTAL_COUNT_OF_TOP_HIT) / $TOTAL_READS_IN_INPUT" | \
       bc -l | awk '{printf "%.3f\n", $0}' | tee '~{reads_basename}.count.~{ref_basename}.pct_total_reads_mapped.txt' )
+
+      if [ $TOTAL_COUNT_OF_TOP_HIT -ne 0 ]; then
+        PCT_TOP_HIT_OF_TOTAL_READS=$( echo "scale=3; 100 * ($TOTAL_COUNT_OF_TOP_HIT / $TOTAL_READS_IN_INPUT)" | \
+          bc -l | awk '{printf "%.3f\n", $0}' | tee '~{reads_basename}.count.~{ref_basename}.pct_top_hit_of_total_reads.txt' )
+      else
+        echo "PCT_TOP_HIT_OF_TOTAL_READS cannot be calculated: there were no mapping hits, or no reads"
+        PCT_TOP_HIT_OF_TOTAL_READS=$( echo "null" | tee '~{reads_basename}.count.~{ref_basename}.pct_top_hit_of_total_reads.txt')
+      fi
     fi
   >>>
 
@@ -493,8 +501,9 @@ task align_and_count {
     Int    reads_mapped_top_hit = read_int("TOTAL_COUNT_OF_TOP_HIT")
     Int    reads_mapped         = read_int("TOTAL_COUNT_OF_LESSER_HITS") + read_int("TOTAL_COUNT_OF_TOP_HIT")
 
-    String pct_total_reads_mapped    = read_string('~{reads_basename}.count.~{ref_basename}.pct_total_reads_mapped.txt')
-    String pct_lesser_hits_of_mapped = read_string('~{reads_basename}.count.~{ref_basename}.pct_lesser_hits_of_mapped.txt')
+    String pct_total_reads_mapped     = read_string('~{reads_basename}.count.~{ref_basename}.pct_total_reads_mapped.txt')
+    String pct_top_hit_of_total_reads = read_string('~{reads_basename}.count.~{ref_basename}.pct_top_hit_of_total_reads.txt')
+    String pct_lesser_hits_of_mapped  = read_string('~{reads_basename}.count.~{ref_basename}.pct_lesser_hits_of_mapped.txt')
     
     String viralngs_version = read_string("VERSION")
   }

diff --git a/pipes/WDL/workflows/align_and_count.wdl b/pipes/WDL/workflows/align_and_count.wdl
@@ -26,8 +26,9 @@ workflow align_and_count_report {
         File   report_top_hits           = align_and_count.report_top_hits
         String tophit                    = align_and_count.top_hit_id
 
-        String pct_mapped_of_total_reads = align_and_count.pct_total_reads_mapped
-        String pct_mapped_to_lesser_hits = align_and_count.pct_lesser_hits_of_mapped
+        String pct_mapped_of_total_reads  = align_and_count.pct_total_reads_mapped
+        String pct_top_hit_of_total_reads = align_and_count.pct_top_hit_of_total_reads
+        String pct_mapped_to_lesser_hits  = align_and_count.pct_lesser_hits_of_mapped
 
         String viral_core_version        = align_and_count.viralngs_version
     }