feat: deduplicate with UMIs

BALSAMIC/assets/scripts/cap_base_quality_in_bam.py

@@ -0,0 +1,37 @@
+import click
+import pysam
+import numpy as np
+
+
+@click.command()
+@click.argument("input_bam", type=click.Path(exists=True))
+@click.argument("output_bam", type=click.Path())
+@click.option(
+    "--max-quality",
+    default=70,
+    type=int,
+    help="Maximum quality value to cap to.",
+)
+def cap_base_qualities(input_bam: str, output_bam: str, max_quality: int):
+    """
+    Cap the base qualities in a BAM file.
+
+    Args:
+        input_bam (str): Input BAM file path.
+        output_bam (str): Output BAM file path.
+        max_quality (int): Maximum quality value to cap to.
+    """
+    # Open input BAM file for reading
+    samfile = pysam.AlignmentFile(input_bam, "rb")
+    out_bam = pysam.AlignmentFile(output_bam, "wb", header=samfile.header)
+    for read in samfile.fetch():
+        qualities = np.array(read.query_qualities)
+        capped_qualities = np.minimum(qualities, max_quality)
+        # Update the base qualities in the read
+        read.query_qualities = capped_qualities.tolist()
+        # Write the modified read to the output BAM file
+        out_bam.write(read)
+
+
+if __name__ == "__main__":
+    cap_base_qualities()

BALSAMIC/assets/scripts/collect_qc_metrics.py

@@ -197,7 +197,6 @@ def get_sample_id(multiqc_key: str) -> str:
     Returns
         str: The extracted sample ID with the ACCXXXXXX format.
     """
-
     if "_align_sort_" in multiqc_key:
         return multiqc_key.split("_")[0]
     return multiqc_key.split(".")[1].split("_")[0]

BALSAMIC/assets/scripts/postprocess_fixmate_bam.awk

@@ -0,0 +1,32 @@
+#!/usr/bin/awk -f
+
+BEGIN { OFS = "\t" }
+/^@/ { print; next }
+{
+    flag = $2
+
+    # If the mate is unmapped, remove the MC tag if it's "*"
+    if (and(flag, 8) != 0) {
+        for (i = 12; i <= NF; i++) {
+            if ($i ~ /^MC:Z:\*$/) {
+                $i = ""
+            }
+        }
+    }
+
+    # Check if any of the specific bitwise flags are set (2, 8, 32, 64, 128)
+    if (and(flag, 2 + 8 + 32 + 64 + 128) != 0) {
+        # Add mate unmapped flag if the mate is unmapped
+        if ($7 == "*" || and(flag, 8) != 0) {
+            flag = or(flag, 8)
+        }
+
+        # Ensure the read paired flag (1) is set if any of these are present
+        flag = or(flag, 1)
+    }
+
+    # Set the modified flag
+    $2 = flag
+
+    print
+}

BALSAMIC/commands/config/case.py

@@ -9,7 +9,6 @@
 
 from BALSAMIC import __version__ as balsamic_version
 from BALSAMIC.commands.options import (
-    OPTION_ADAPTER_TRIM,
     OPTION_ANALYSIS_DIR,
     OPTION_ANALYSIS_WORKFLOW,
     OPTION_BACKGROUND_VARIANTS,
@@ -32,15 +31,12 @@
     OPTION_NORMAL_SAMPLE_NAME,
     OPTION_PANEL_BED,
     OPTION_PON_CNN,
-    OPTION_QUALITY_TRIM,
     OPTION_SENTIEON_INSTALL_DIR,
     OPTION_SENTIEON_LICENSE,
     OPTION_SWEGEN_SNV,
     OPTION_SWEGEN_SV,
     OPTION_TUMOR_SAMPLE_NAME,
-    OPTION_UMI,
     OPTION_UMI_MIN_READS,
-    OPTION_UMI_TRIM_LENGTH,
 )
 from BALSAMIC.constants.analysis import BIOINFO_TOOL_ENV, AnalysisWorkflow, Gender
 from BALSAMIC.constants.cache import GenomeVersion
@@ -66,7 +62,6 @@
 
 
 @click.command("case", short_help="Create a sample config file from input sample data")
-@OPTION_ADAPTER_TRIM
 @OPTION_ANALYSIS_DIR
 @OPTION_ANALYSIS_WORKFLOW
 @OPTION_BACKGROUND_VARIANTS
@@ -89,19 +84,15 @@
 @OPTION_NORMAL_SAMPLE_NAME
 @OPTION_PANEL_BED
 @OPTION_PON_CNN
-@OPTION_QUALITY_TRIM
 @OPTION_SENTIEON_INSTALL_DIR
 @OPTION_SENTIEON_LICENSE
 @OPTION_SWEGEN_SNV
 @OPTION_SWEGEN_SV
 @OPTION_TUMOR_SAMPLE_NAME
-@OPTION_UMI
-@OPTION_UMI_TRIM_LENGTH
 @OPTION_UMI_MIN_READS
 @click.pass_context
 def case_config(
     context: click.Context,
-    adapter_trim: bool,
     analysis_dir: Path,
     analysis_workflow: AnalysisWorkflow,
     background_variants: Path,
@@ -124,14 +115,11 @@ def case_config(
     normal_sample_name: str,
     panel_bed: Path,
     pon_cnn: Path,
-    quality_trim: bool,
     sentieon_install_dir: Path,
     sentieon_license: str,
     swegen_snv: Path,
     swegen_sv: Path,
     tumor_sample_name: str,
-    umi: bool,
-    umi_trim_length: int,
     umi_min_reads: str | None,
 ):
     references_path: Path = Path(balsamic_cache, cache_version, genome_version)
@@ -205,12 +193,6 @@ def case_config(
             "dnascope_model": SENTIEON_DNASCOPE_MODEL.as_posix(),
             "tnscope_model": SENTIEON_TNSCOPE_MODEL.as_posix(),
         },
-        QC={
-            "quality_trim": quality_trim,
-            "adapter_trim": adapter_trim,
-            "umi_trim": umi if panel_bed else False,
-            "umi_trim_length": umi_trim_length,
-        },
         analysis={
             "case_id": case_id,
             "gender": gender,

BALSAMIC/commands/config/pon.py

@@ -8,7 +8,6 @@
 
 from BALSAMIC import __version__ as balsamic_version
 from BALSAMIC.commands.options import (
-    OPTION_ADAPTER_TRIM,
     OPTION_ANALYSIS_DIR,
     OPTION_BALSAMIC_CACHE,
     OPTION_CACHE_VERSION,
@@ -21,9 +20,6 @@
     OPTION_PANEL_BED,
     OPTION_PON_VERSION,
     OPTION_PON_WORKFLOW,
-    OPTION_QUALITY_TRIM,
-    OPTION_UMI,
-    OPTION_UMI_TRIM_LENGTH,
 )
 from BALSAMIC.constants.analysis import BIOINFO_TOOL_ENV, PONWorkflow
 from BALSAMIC.constants.cache import GenomeVersion
@@ -47,7 +43,6 @@
 
 
 @click.command("pon", short_help="Create a sample config file for PON analysis")
-@OPTION_ADAPTER_TRIM
 @OPTION_ANALYSIS_DIR
 @OPTION_BALSAMIC_CACHE
 @OPTION_CACHE_VERSION
@@ -60,13 +55,9 @@
 @OPTION_PANEL_BED
 @OPTION_PON_WORKFLOW
 @OPTION_PON_VERSION
-@OPTION_QUALITY_TRIM
-@OPTION_UMI
-@OPTION_UMI_TRIM_LENGTH
 @click.pass_context
 def pon_config(
     context: click.Context,
-    adapter_trim: bool,
     analysis_dir: Path,
     balsamic_cache: Path,
     cache_version: str,
@@ -78,9 +69,6 @@ def pon_config(
     sentieon_license: str,
     panel_bed: Path,
     pon_workflow: PONWorkflow,
-    quality_trim: bool,
-    umi: bool,
-    umi_trim_length: bool,
     version: str,
 ):
     references_path: Path = Path(balsamic_cache, cache_version, genome_version)
@@ -122,12 +110,6 @@ def pon_config(
             "dnascope_model": SENTIEON_DNASCOPE_MODEL.as_posix(),
             "tnscope_model": SENTIEON_TNSCOPE_MODEL.as_posix(),
         },
-        QC={
-            "adapter_trim": adapter_trim,
-            "quality_trim": quality_trim,
-            "umi_trim": umi if panel_bed else False,
-            "umi_trim_length": umi_trim_length,
-        },
         analysis={
             "case_id": case_id,
             "analysis_dir": analysis_dir,

BALSAMIC/commands/options.py

@@ -29,14 +29,6 @@
     validate_umi_min_reads,
 )
 
-OPTION_ADAPTER_TRIM = click.option(
-    "--adapter-trim/--no-adapter-trim",
-    default=True,
-    show_default=True,
-    is_flag=True,
-    help="Trim adapters from reads in FASTQ file",
-)
-
 OPTION_ANALYSIS_DIR = click.option(
     "--analysis-dir",
     type=click.Path(exists=True, resolve_path=True),
@@ -319,14 +311,6 @@
     help="Print list of analysis files. Otherwise only final count will be printed.",
 )
 
-OPTION_QUALITY_TRIM = click.option(
-    "--quality-trim/--no-quality-trim",
-    default=True,
-    show_default=True,
-    is_flag=True,
-    help="Trim low quality reads in FASTQ file",
-)
-
 OPTION_QUIET = click.option(
     "-q",
     "--quiet",
@@ -427,22 +411,6 @@
     help="Tumor sample name",
 )
 
-OPTION_UMI = click.option(
-    "--umi/--no-umi",
-    default=True,
-    show_default=True,
-    is_flag=True,
-    help="UMI processing steps for samples with UMI tags. For WGS cases, UMI is always disabled.",
-)
-
-OPTION_UMI_TRIM_LENGTH = click.option(
-    "--umi-trim-length",
-    default=5,
-    show_default=True,
-    type=click.INT,
-    help="Trim N bases from reads in FASTQ file",
-)
-
 OPTION_UMI_MIN_READS = click.option(
     "--umi-min-reads",
     type=click.STRING,

BALSAMIC/constants/analysis.py

@@ -124,6 +124,7 @@ class BioinfoTools(StrEnum):
     VCFANNO: str = "vcfanno"
     CADD: str = "cadd"
     PURECN: str = "purecn"
+    PYSAM: str = "pysam"
 
 
 class FastqName(StrEnum):
@@ -177,4 +178,5 @@ class PonParams:
     BioinfoTools.CADD: DockerContainers.CADD,
     BioinfoTools.PURECN: DockerContainers.PURECN,
     BioinfoTools.GATK: DockerContainers.GATK,
+    BioinfoTools.PYSAM: DockerContainers.PYTHON_3,
 }

BALSAMIC/constants/cluster_analysis.json

@@ -19,6 +19,10 @@
 		"time": "03:00:00",
 		"n": 12
 	},
+	"cap_base_quality": {
+		"time": "05:00:00",
+		"n": 10
+	},
 	"finalize_gens_outputfiles": {
 		"time": "01:00:00",
 		"n": 2
@@ -132,10 +136,18 @@
 		"time": "02:30:00",
 		"n": 8
 	},
-	"samtools_sort_index": {
-		"time": "02:30:00",
+	"samtools_fixmate": {
+		"time": "04:30:00",
 		"n": 16
 	},
+	"sentieon_qc_metrics": {
+		"time": "06:00:00",
+		"n": 8
+	},
+	"sentieon_plot_qc_metrics": {
+		"time": "03:00:00",
+		"n": 2
+	},
 	"sentieon_DNAscope": {
 		"time": "24:00:00",
 		"n": 36