Add FASTQ header cleaning as the first pipeline step

.github/workflows/nf-test.yml

@@ -36,6 +36,7 @@ jobs:
           - "only_trimming/skip_trimming"
           - "only_alignment/alignment"
           - "only_moveumitoheader"
+          - "only_cleanheaders"
     steps:
       - name: Check out pipeline code
         uses: actions/checkout@v3

Documentation.md

@@ -37,6 +37,12 @@ MultiQC report of important quality control metrics. Produced by MULTIQC process
 
 ## Commonly used parameters
 
+**Cleaning spaces from fastq headers**
+
+Some fastq files carry spaces in their header lines (for example `@READID 1:N:0:INDEX`). Tools further down the pipeline truncate read names at the first space, which can break UMI extraction and read-name matching after alignment.
+
+To enable this option set `clean_fastq_headers = true`. Every space in a header line is replaced with an underscore, and the cleaned reads are used for all subsequent steps. This runs before any other analysis, including UMI moving and trimming.
+
 **Moving UMI from fastq reads to read header**
 
 If you are analysing demultiplexed sample files, then depending on where you have sourced your fastq from, the UMI and experimental barcode might still be present at the 5' end of reads, which will cause errors in downstream analysis. Note this is common for historical public data downloaded from ArrayExpress, for example.

conf/logic.config

@@ -7,6 +7,7 @@
 
 // Defaults
 params {
+    run_clean_fastq_headers = false
     run_genome_prep         = true
     run_input_check         = true
     run_move_umi_to_header  = false
@@ -20,6 +21,7 @@ params {
 }
 
 // Set other logic
+if(params.clean_fastq_headers) { params.run_clean_fastq_headers = true }
 if(params.move_umi_to_header) { params.run_move_umi_to_header = true }
 if(params.skip_umi_dedupe)    { params.run_umi_dedup = false         }
 
@@ -36,6 +38,7 @@ if(params.only_input) {
 }
 
 if(params.only_genome) {
+    params.run_clean_fastq_headers = false
     params.run_input_check        = false
     params.run_trim_galore_fastqc = false
     params.run_alignment          = false

conf/modules.config

@@ -196,6 +196,18 @@ if(params.run_genome_prep) {
     FASTQ + TRIMMING
 ========================================================================================
 */
+if(params.run_clean_fastq_headers) {
+    process {
+        withName: 'CLIPSEQ:CLEAN_FASTQ_HEADERS' {
+            publishDir = [
+                path: { "${params.outdir}/01_prealign/clean_fastq_headers" },
+                mode: "${params.publish_dir_mode}",
+                enabled: false
+            ]
+        }
+    }
+}
+
 if(params.run_move_umi_to_header) {
     process {
         withName: 'CLIPSEQ:UMITOOLS_EXTRACT' {

main.nf

@@ -94,6 +94,7 @@ ch_multiqc_config = file("$projectDir/assets/multiqc_config.yml", checkIfExists:
 // MODULEs
 //
 
+include { CLEAN_FASTQ_HEADERS } from './modules/local/clean_fastq_headers/main'
 include { MULTIQC } from './modules/local/multiqc'
 include { GET_CROSSLINKS as CALC_SMRNA_K1_CROSSLINKS      } from './modules/local/get_crosslinks'
 include { GET_CROSSLINKS as CALC_GENOME_CROSSLINKS        } from './modules/local/get_crosslinks'
@@ -279,6 +280,17 @@ workflow CLIPSEQ {
     }
     //EXAMPLE CHANNEL STRUCT: [[id:h3k27me3_R1, group:h3k27me3, replicate:1, single_end:false], [FASTQ]]
     //ch_fastq | view
+    if(params.run_clean_fastq_headers) {
+        /*
+        * MODULE: Replace spaces in FASTQ headers with underscores before any analysis
+        */
+        CLEAN_FASTQ_HEADERS (
+            ch_fastq
+        )
+        ch_versions = ch_versions.mix(CLEAN_FASTQ_HEADERS.out.versions)
+        ch_fastq    = CLEAN_FASTQ_HEADERS.out.reads
+    }
+
     if(params.encode_eclip){
         ENCODE_MOVEUMI (
             ch_fastq

modules/local/clean_fastq_headers/main.nf

@@ -0,0 +1,22 @@
+process CLEAN_FASTQ_HEADERS {
+    tag "$meta.id"
+    label "process_single"
+
+    conda "bioconda::biopython=1.78 pigz=2.6"
+    container "quay.io/biocontainers/mulled-v2-877c4e5a8fad685ea5bde487e04924ac447923b9:b7daa641364165419b9a87d9988bc803f913c5b6-0"
+
+    input:
+    tuple val(meta), path(reads)
+
+    output:
+    tuple val(meta), path("*.clean.fastq.gz"), emit: reads
+    path  "versions.yml"                     , emit: versions
+
+    when:
+    task.ext.when == null || task.ext.when
+
+    shell:
+    prefix       = task.ext.prefix ?: "${meta.id}"
+    process_name = task.process
+    template 'clean_fastq_headers.py'
+}

modules/local/clean_fastq_headers/meta.yml

@@ -0,0 +1,37 @@
+name: clean_fastq_headers
+description: Replace spaces in FASTQ header lines with underscores before any analysis
+tools:
+  - pigz:
+      description: |
+        A parallel implementation of gzip for modern multi-processor,
+        multi-core machines.
+      homepage: https://zlib.net/pigz/
+      documentation: https://zlib.net/pigz/pigz.pdf
+      licence: ["Zlib"]
+input:
+  - meta:
+      type: map
+      description: |
+        Groovy Map containing sample information
+        e.g. [ id:'test', single_end:false ]
+  - reads:
+      type: file
+      description: A single FASTQ file (this pipeline is single-end only)
+      pattern: "*.{fastq,fastq.gz}"
+
+output:
+  - meta:
+      type: map
+      description: |
+        Groovy Map containing sample information
+        e.g. [ id:'test', single_end:false ]
+  - reads:
+      type: file
+      description: FASTQ file with spaces in the header lines replaced by underscores
+      pattern: "*.clean.fastq.gz"
+  - versions:
+      type: file
+      description: File containing software versions
+      pattern: "versions.yml"
+authors:
+  - "@Chromojones"

modules/local/clean_fastq_headers/templates/clean_fastq_headers.py

@@ -0,0 +1,76 @@
+#!/usr/bin/env python3
+
+import platform
+import subprocess
+
+
+def process_fastq(input_file, output_file, threads):
+    """Process FASTQ file to remove spaces from headers using pigz."""
+    is_gzipped = input_file.endswith('.gz')
+
+    in_file = None
+    if is_gzipped:
+        decompress = subprocess.Popen(['pigz', '-p', str(threads), '-dc', input_file],
+                                      stdout=subprocess.PIPE)
+    else:
+        in_file = open(input_file, 'rb')
+        decompress = subprocess.Popen(['cat'], stdin=in_file, stdout=subprocess.PIPE)
+
+    with open(output_file, 'wb') as out_handle:
+        compress = subprocess.Popen(['pigz', '-p', str(threads), '-c'],
+                                    stdin=subprocess.PIPE, stdout=out_handle)
+        try:
+            line_count = 0
+            for line in decompress.stdout:
+                if line_count % 4 == 0:  # Header line
+                    line = line.decode().strip().replace(' ', '_').encode() + b'\n'
+                compress.stdin.write(line)
+                line_count += 1
+        finally:
+            decompress.stdout.close()
+            decompress.wait()
+            if in_file is not None:
+                in_file.close()
+            compress.stdin.close()
+            compress.wait()
+
+    # Both ends of the pipe must be checked. pigz reports a corrupt or
+    # truncated input by exiting non-zero after emitting a partial stream, so
+    # without this the task "succeeds" and hands a silently truncated FASTQ to
+    # the rest of the pipeline.
+    if decompress.returncode != 0:
+        raise RuntimeError(
+            "reading %s failed (exit %d) - the file is likely truncated or corrupt"
+            % (input_file, decompress.returncode))
+    if compress.returncode != 0:
+        raise RuntimeError(
+            "writing %s failed (exit %d)" % (output_file, compress.returncode))
+
+    if line_count % 4 != 0:
+        raise RuntimeError(
+            "%s ended mid-record after %d lines - the file is likely truncated"
+            % (input_file, line_count))
+
+
+def pigz_version():
+    """pigz reports its version on stderr on some builds, stdout on others."""
+    proc = subprocess.run(['pigz', '--version'], stdout=subprocess.PIPE,
+                          stderr=subprocess.STDOUT)
+    return proc.stdout.decode().strip().replace('pigz ', '')
+
+
+# This pipeline is single-end only: single_end is parsed into meta but never
+# read downstream, and each sample is one FASTQ. Fail loudly rather than
+# quietly emit files the rest of the pipeline cannot use.
+reads = "!{reads}".replace('[', '').replace(']', '').replace(',', ' ').split()
+if len(reads) != 1:
+    raise SystemExit(
+        "CLEAN_FASTQ_HEADERS expects exactly one FASTQ per sample, got %d (%s)"
+        % (len(reads), ", ".join(reads)))
+
+process_fastq(reads[0], "!{prefix}.clean.fastq.gz", !{task.cpus})
+
+with open("versions.yml", "w") as out_f:
+    out_f.write("!{process_name}" + ":\n")
+    out_f.write("    python: " + platform.python_version() + "\n")
+    out_f.write("    pigz: " + pigz_version() + "\n")

nextflow.config

@@ -14,7 +14,7 @@ params {
     publish_dir_mode = 'symlink'
     monochrome_logs  = false
     debug            = false
-    ignore_params    = "run_genome_prep,run_input_check,run_trim_galore_fastqc,run_alignment,run_read_filter,run_umi_dedup,run_calc_crosslinks,run_peak_calling,run_reporting,only_input,only_genome,only_trimming,only_alignment,only_filtering,only_dedup,only_crosslinks,only_peakcalling"
+    ignore_params    = "run_clean_fastq_headers,run_genome_prep,run_input_check,run_trim_galore_fastqc,run_alignment,run_read_filter,run_umi_dedup,run_calc_crosslinks,run_peak_calling,run_reporting,only_input,only_genome,only_trimming,only_alignment,only_filtering,only_dedup,only_crosslinks,only_peakcalling"
     multiqc_title    = null
 
     // Max resource options
@@ -70,6 +70,7 @@ params {
 
     // Pipeline params
     crosslink_position  = "start"
+    clean_fastq_headers = false
     encode_eclip        = false
     move_umi_to_header  = false
     umi_header_format   = null

schema/clipseq.json

@@ -303,6 +303,11 @@
                     "description": "Format of the UMI header to extract.",
                     "type": "string"
                 },
+                "clean_fastq_headers": {
+                    "name": "Clean FASTQ Headers",
+                    "description": "Replace spaces in FASTQ header lines with underscores before any analysis.",
+                    "type": "boolean"
+                },
                 "skip_umi_dedupe": {
                     "name": "Skip UMI Deduplication",
                     "description": "Skip the UMI deduplication of BAMs.",

tests/nf_test/only_cleanheaders/main.nf.test

@@ -0,0 +1,28 @@
+nextflow_pipeline {
+
+    name "only_cleanheaders"
+    script "main.nf"
+
+    test("run_test") {
+        when {
+            params {
+                outdir = "$outputDir"
+                only_trimming = true
+                save_trimmed = true
+                clean_fastq_headers = true
+            }
+        }
+
+        then {
+            assert workflow.success
+
+            // The cleaned fastqs are intermediates and are not published, so
+            // assert the process ran. TrimGalore re-links its input to
+            // "${prefix}.fastq.gz", so downstream file names are unchanged.
+            assert workflow.trace.succeeded().any { it.process ==~ /.*CLEAN_FASTQ_HEADERS.*/ }
+
+            assert new File("$outputDir/01_prealign/trimgalore/sample_1_R1.fastq.gz_trimming_report.txt").exists()
+            assert new File("$outputDir/01_prealign/trimgalore/sample_1_R1_trimmed.fq.gz").exists()
+        }
+    }
+}