feat: report mean gap-compressed identity for primary and supplementary alignments

williamrowell · williamrowell · commit fd09bc969028 · 2025-09-04T09:01:03.000-07:00
docs: Update documentation for output of bam_statistics task.

feat: add read length n50 calculation

Rebase submodule.
diff --git a/docs/bam_statistics.md b/docs/bam_statistics.md
@@ -0,0 +1,37 @@
+# bam_stats outputs
+
+## `bam_statistics`
+
+A compressed TSV file with a row for each record in the haplotagged BAM and the following columns:
+
+- movie name
+- read name
+- read length
+- Phred scaled read quality
+- alignment type (unmapped, primary, supplementary; because supplementary alignments are included, reads may appear on multiple rows)
+- mapping quality (`MAPQ`), if mapped
+- gap-compressed identity (`mg`), if mapped
+
+## `read_length_plot`
+
+A histogram of read lengths, using only records marked `prim` or `unmapped`.
+
+## `read_quality_plot`
+
+A histogram of read qualities, using only records marked `prim` or `unmapped`. This output is only generated if the input BAMs contain the `rq` tag.
+
+## `mapq_distribution_plot`, `mg_distribution_plot`
+
+A histogram of mapping qualities and gap-compressed identities, respectively.
+
+## `stat_num_reads`, `stat_read_length_mean`, `stat_read_length_median`, `stat_read_length_n50`, `stat_read_quality_mean`, `stat_read_quality_median`
+
+Statistics computed using only records marked `prim` or `unmapped`.
+
+## `stat_mapped_read_count`, `stat_mapped_percent`
+
+Count of primary alignments, and primary alignments as a percentage of total reads.
+
+## `stat_mean_gap_compressed_identity`
+
+Mean gap-compressed identity of primary and supplementary alignments.
diff --git a/docs/bam_stats.md b/docs/bam_stats.md
diff --git a/docs/family.md b/docs/family.md
@@ -169,10 +169,12 @@ The `Sample` struct contains sample specific data and metadata. The struct has t
 | Array\[String\] | stat_num_reads | Number of reads |  |
 | Array\[String\] | stat_read_length_mean | Mean read length |  |
 | Array\[String\] | stat_read_length_median | Median read length |  |
+| Array\[String\] | stat_read_length_n50 | Read length N50 |  |
 | Array\[String\] | stat_read_quality_mean | Mean read quality |  |
 | Array\[String\] | stat_read_quality_median | Median read quality |  |
 | Array\[String\] | stat_mapped_read_count | Count of reads mapped to reference |  |
 | Array\[String\] | stat_mapped_percent | Percent of reads mapped to reference |  |
+| Array\[String\] | stat_mean_gap_compressed_identity | Mean gap-compressed identity |  |
 | Array\[String\] | inferred_sex | Inferred sex | Sex is inferred based on relative depth of chrY alignments. |
 | Array\[String\] | stat_mean_depth | Mean depth | |
 
diff --git a/docs/singleton.md b/docs/singleton.md
@@ -125,10 +125,12 @@ flowchart TD
 | String | stat_num_reads | Number of reads |  |
 | String | stat_read_length_mean | Mean read length |  |
 | String | stat_read_length_median | Median read length |  |
+| String | stat_read_length_n50 | Read length N50 |  |
 | String | stat_read_quality_mean | Mean read quality |  |
 | String | stat_read_quality_median | Median read quality |  |
 | String | stat_mapped_read_count | Count of reads mapped to reference |  |
 | String | stat_mapped_percent | Percent of reads mapped to reference |  |
+| String | stat_mean_gap_compressed_identity | Mean gap-compressed identity |  |
 | String | inferred_sex | Inferred sex | Sex is inferred based on relative depth of chrY alignments. |
 | String | stat_mean_depth | Mean depth | |
 
diff --git a/wdl-ci.config.json b/wdl-ci.config.json
@@ -269,7 +269,7 @@
       "tasks": {
         "bam_stats": {
           "key": "bam_stats",
-          "digest": "3gn3hurjmdifhucjdrnykji4w4cf7yjq",
+          "digest": "n7yzgzllk24zqced7wsgijsmdjoot66t",
           "tests": [
             {
               "inputs": {
@@ -339,6 +339,12 @@
                     "compare_string"
                   ]
                 },
+                "stat_read_length_n50": {
+                  "value": "16945",
+                  "test_tasks": [
+                    "compare_string"
+                  ]
+                },
                 "stat_read_quality_mean": {
                   "value": "35.91",
                   "test_tasks": [
@@ -362,6 +368,12 @@
                   "test_tasks": [
                     "compare_string"
                   ]
+                },
+                "stat_mean_gap_compressed_identity": {
+                  "value": "99.77",
+                  "test_tasks": [
+                    "compare_string"
+                  ]
                 }
               }
             }
diff --git a/workflows/downstream/downstream.wdl b/workflows/downstream/downstream.wdl
@@ -212,19 +212,21 @@ workflow downstream {
     String stat_phase_block_ng50          = hiphase.stat_phase_block_ng50
 
     # bam stats
-    File   bam_statistics           = bam_stats.bam_statistics
-    File   read_length_plot         = bam_stats.read_length_plot
-    File?  read_quality_plot        = bam_stats.read_quality_plot
-    File   mapq_distribution_plot   = bam_stats.mapq_distribution_plot
-    File   mg_distribution_plot     = bam_stats.mg_distribution_plot
-    String stat_num_reads           = bam_stats.stat_num_reads
-    String stat_read_length_mean    = bam_stats.stat_read_length_mean
-    String stat_read_length_median  = bam_stats.stat_read_length_median
-    String stat_read_quality_mean   = bam_stats.stat_read_quality_mean
-    String stat_read_quality_median = bam_stats.stat_read_quality_median
-    String stat_mapped_read_count   = bam_stats.stat_mapped_read_count
-    String stat_mapped_percent      = bam_stats.stat_mapped_percent
-    File   trgt_coverage_dropouts   = coverage_dropouts.dropouts
+    File   bam_statistics                    = bam_stats.bam_statistics
+    File   read_length_plot                  = bam_stats.read_length_plot
+    File?  read_quality_plot                 = bam_stats.read_quality_plot
+    File   mapq_distribution_plot            = bam_stats.mapq_distribution_plot
+    File   mg_distribution_plot              = bam_stats.mg_distribution_plot
+    String stat_num_reads                    = bam_stats.stat_num_reads
+    String stat_read_length_mean             = bam_stats.stat_read_length_mean
+    String stat_read_length_median           = bam_stats.stat_read_length_median
+    String stat_read_length_n50              = bam_stats.stat_read_length_n50
+    String stat_read_quality_mean            = bam_stats.stat_read_quality_mean
+    String stat_read_quality_median          = bam_stats.stat_read_quality_median
+    String stat_mapped_read_count            = bam_stats.stat_mapped_read_count
+    String stat_mapped_percent               = bam_stats.stat_mapped_percent
+    String stat_mean_gap_compressed_identity = bam_stats.stat_mean_gap_compressed_identity
+    File   trgt_coverage_dropouts            = coverage_dropouts.dropouts
 
     # small variant stats
     File   small_variant_stats     = bcftools_stats_roh_small_variants.stats
diff --git a/workflows/family.wdl b/workflows/family.wdl
@@ -224,10 +224,12 @@ workflow humanwgs_family {
     'num_reads': downstream.stat_num_reads,
     'read_length_mean': downstream.stat_read_length_mean,
     'read_length_median': downstream.stat_read_length_median,
+    'read_length_n50': downstream.stat_read_length_n50,
     'read_quality_mean': downstream.stat_read_quality_mean,
     'read_quality_median': downstream.stat_read_quality_median,
     'mapped_read_count': downstream.stat_mapped_read_count,
     'mapped_percent': downstream.stat_mapped_percent,
+    'mean_gap_compressed_identity': downstream.stat_mean_gap_compressed_identity,
     'mean_depth': upstream.stat_mean_depth,
     'inferred_sex': upstream.inferred_sex,
     'stat_phased_basepairs': downstream.stat_phased_basepairs,
@@ -267,18 +269,20 @@ workflow humanwgs_family {
     File  msg_file           = consolidate_stats.messages
 
     # bam stats
-    Array[File]   bam_statistics           = downstream.bam_statistics
-    Array[File]   read_length_plot         = downstream.read_length_plot
-    Array[File?]  read_quality_plot        = downstream.read_quality_plot
-    Array[File]   mapq_distribution_plot   = downstream.mapq_distribution_plot
-    Array[File]   mg_distribution_plot     = downstream.mg_distribution_plot
-    Array[String] stat_num_reads           = downstream.stat_num_reads
-    Array[String] stat_read_length_mean    = downstream.stat_read_length_mean
-    Array[String] stat_read_length_median  = downstream.stat_read_length_median
-    Array[String] stat_read_quality_mean   = downstream.stat_read_quality_mean
-    Array[String] stat_read_quality_median = downstream.stat_read_quality_median
-    Array[String] stat_mapped_read_count   = downstream.stat_mapped_read_count
-    Array[String] stat_mapped_percent      = downstream.stat_mapped_percent
+    Array[File]   bam_statistics                    = downstream.bam_statistics
+    Array[File]   read_length_plot                  = downstream.read_length_plot
+    Array[File?]  read_quality_plot                 = downstream.read_quality_plot
+    Array[File]   mapq_distribution_plot            = downstream.mapq_distribution_plot
+    Array[File]   mg_distribution_plot              = downstream.mg_distribution_plot
+    Array[String] stat_num_reads                    = downstream.stat_num_reads
+    Array[String] stat_read_length_mean             = downstream.stat_read_length_mean
+    Array[String] stat_read_length_median           = downstream.stat_read_length_median
+    Array[String] stat_read_length_n50              = downstream.stat_read_length_n50
+    Array[String] stat_read_quality_mean            = downstream.stat_read_quality_mean
+    Array[String] stat_read_quality_median          = downstream.stat_read_quality_median
+    Array[String] stat_mapped_read_count            = downstream.stat_mapped_read_count
+    Array[String] stat_mapped_percent               = downstream.stat_mapped_percent
+    Array[String] stat_mean_gap_compressed_identity = downstream.stat_mean_gap_compressed_identity
 
     # merged, haplotagged alignments
     Array[File]   merged_haplotagged_bam       = downstream.merged_haplotagged_bam
diff --git a/workflows/singleton.wdl b/workflows/singleton.wdl
@@ -169,10 +169,12 @@ workflow humanwgs_singleton {
     'num_reads': [downstream.stat_num_reads],
     'read_length_mean': [downstream.stat_read_length_mean],
     'read_length_median': [downstream.stat_read_length_median],
+    'read_length_n50': [downstream.stat_read_length_n50],
     'read_quality_mean': [downstream.stat_read_quality_mean],
     'read_quality_median': [downstream.stat_read_quality_median],
     'mapped_read_count': [downstream.stat_mapped_read_count],
     'mapped_percent': [downstream.stat_mapped_percent],
+    'mean_gap_compressed_identity': [downstream.stat_mean_gap_compressed_identity],
     'mean_depth': [upstream.stat_mean_depth],
     'inferred_sex': [upstream.inferred_sex],
     'stat_phased_basepairs': [downstream.stat_phased_basepairs],
@@ -211,18 +213,20 @@ workflow humanwgs_singleton {
     File msg_file   = consolidate_stats.messages
 
     # bam stats
-    File   bam_statistics           = downstream.bam_statistics
-    File   read_length_plot         = downstream.read_length_plot
-    File?  read_quality_plot        = downstream.read_quality_plot
-    File   mapq_distribution_plot   = downstream.mapq_distribution_plot
-    File   mg_distribution_plot     = downstream.mg_distribution_plot
-    String stat_num_reads           = downstream.stat_num_reads
-    String stat_read_length_mean    = downstream.stat_read_length_mean
-    String stat_read_length_median  = downstream.stat_read_length_median
-    String stat_read_quality_mean   = downstream.stat_read_quality_mean
-    String stat_read_quality_median = downstream.stat_read_quality_median
-    String stat_mapped_read_count   = downstream.stat_mapped_read_count
-    String stat_mapped_percent      = downstream.stat_mapped_percent
+    File   bam_statistics                    = downstream.bam_statistics
+    File   read_length_plot                  = downstream.read_length_plot
+    File?  read_quality_plot                 = downstream.read_quality_plot
+    File   mapq_distribution_plot            = downstream.mapq_distribution_plot
+    File   mg_distribution_plot              = downstream.mg_distribution_plot
+    String stat_num_reads                    = downstream.stat_num_reads
+    String stat_read_length_mean             = downstream.stat_read_length_mean
+    String stat_read_length_median           = downstream.stat_read_length_median
+    String stat_read_length_n50              = downstream.stat_read_length_n50
+    String stat_read_quality_mean            = downstream.stat_read_quality_mean
+    String stat_read_quality_median          = downstream.stat_read_quality_median
+    String stat_mapped_read_count            = downstream.stat_mapped_read_count
+    String stat_mapped_percent               = downstream.stat_mapped_percent
+    String stat_mean_gap_compressed_identity = downstream.stat_mean_gap_compressed_identity
 
     # merged, haplotagged alignments
     File   merged_haplotagged_bam       = downstream.merged_haplotagged_bam
diff --git a/workflows/wdl-common b/workflows/wdl-common
@@ -1 +1 @@
-Subproject commit 4c516f7627a5b9ed3e7b5ea4d4cec5b49958a1f2
+Subproject commit 59f9f5c6b4001c277a176331399992bf9a89ae77

Original file line number	Diff line number	Diff line change
`@@ -269,7 +269,7 @@`
`269`	`269`	`"tasks": {`
`270`	`270`	`"bam_stats": {`
`271`	`271`	`"key": "bam_stats",`
`272`		`- "digest": "3gn3hurjmdifhucjdrnykji4w4cf7yjq",`
	`272`	`+ "digest": "n7yzgzllk24zqced7wsgijsmdjoot66t",`
`273`	`273`	`"tests": [`
`274`	`274`	`{`
`275`	`275`	`"inputs": {`
`@@ -339,6 +339,12 @@`
`339`	`339`	`"compare_string"`
`340`	`340`	`]`
`341`	`341`	`},`
	`342`	`+ "stat_read_length_n50": {`
	`343`	`+ "value": "16945",`
	`344`	`+ "test_tasks": [`
	`345`	`+ "compare_string"`
	`346`	`+ ]`
	`347`	`+ },`
`342`	`348`	`"stat_read_quality_mean": {`
`343`	`349`	`"value": "35.91",`
`344`	`350`	`"test_tasks": [`
`@@ -362,6 +368,12 @@`
`362`	`368`	`"test_tasks": [`
`363`	`369`	`"compare_string"`
`364`	`370`	`]`
	`371`	`+ },`
	`372`	`+ "stat_mean_gap_compressed_identity": {`
	`373`	`+ "value": "99.77",`
	`374`	`+ "test_tasks": [`
	`375`	`+ "compare_string"`
	`376`	`+ ]`
`365`	`377`	`}`
`366`	`378`	`}`
`367`	`379`	`}`