Merge pull request #151 from Joon-Klaps/custom-vcf

Joon-Klaps · web-flow · commit 1f1c6f302f3f · 2024-12-03T10:15:41.000+01:00
Add new module to inculde custom mpileup-vcf file for intra-host analyses
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -18,6 +18,8 @@ Initial release of Joon-Klaps/viralgenie, created with the [nf-core](https://nf-
 - Adding mash-screen output to result table ([#140](https://github.com/Joon-Klaps/viralgenie/pull/140))
 - Add logic to allow samples with no reference hits to be analysed ([#141](https://github.com/Joon-Klaps/viralgenie/pull/141))
 - Add visualisation for hybrid scaffold ([#143](https://github.com/Joon-Klaps/viralgenie/pull/143))
+- Add new module to inculde custom mpileup-vcf file for intra-host analyses ([#151](https://github.com/Joon-Klaps/viralgenie/pull/151))
+
 
 ### `Fixed`
 
diff --git a/bin/custom_mpileup.py b/bin/custom_mpileup.py
@@ -0,0 +1,110 @@
+#!/usr/bin/env python
+
+"""Provide a command line tool to filter generate a custom tsv -csv file from a bam"""
+
+import pysam
+import pysamstats
+import csv
+import numpy as np
+from numpy.typing import NDArray
+import argparse
+from pathlib import Path
+import logging
+
+# Initialize logger
+logger = logging.getLogger(__name__)
+logging.basicConfig(level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s")
+
+
+def parse_args(argv=None):
+    """Define and immediately parse command line arguments."""
+    parser = argparse.ArgumentParser(
+        description="Provide a command line tool to create a custom summary of mpileup results from a bam/cram/sam file",
+        epilog="Example: python custom_mpileup.py --clusters_summary file1,file2,file3,... ",
+    )
+    parser = argparse.ArgumentParser(description="Custom mpileup processing script.")
+    parser.add_argument("--alignment", type=Path, help="Input BAM file prefix")
+    parser.add_argument("--reference", type=Path, help="Reference FASTA file")
+    parser.add_argument("--prefix", type=str, help="Name of the output file")
+    parser.add_argument(
+        "-l",
+        "--log-level",
+        help="The desired log level (default WARNING).",
+        choices=("CRITICAL", "ERROR", "WARNING", "INFO", "DEBUG"),
+        default="INFO",
+    )
+    return parser.parse_args(argv)
+
+
+def process_mpileup(filename: Path, reference: Path) -> NDArray:
+    """
+    Process mpileup data using numpy vectorized operations.
+
+    Args:
+        filename: Path to the alignment file (BAM/CRAM/SAM)
+        reference: Path to the reference FASTA file
+
+    Returns:
+        NDArray: Array with columns [position, A, C, G, T, insertions, deletions, consensus]
+    """
+    # Initialize FASTA file properly
+    fasta = pysam.FastaFile(str(reference))
+
+    alignment_file = pysam.AlignmentFile(filename, "rc" if filename.suffix == ".cram" else "rb", reference_filename=str(reference))
+
+    # Convert generator to structured numpy array
+    data = np.array(
+        [
+            (r["pos"], r["ref"], r["A"], r["C"], r["G"], r["T"], r["insertions"], r["deletions"], "N")
+            for r in pysamstats.stat_variation(alignment_file, fafile=fasta)
+        ],
+        dtype=[
+            ("pos", int),
+            ("ref", "U1"),
+            ("A", int),
+            ("C", int),
+            ("G", int),
+            ("T", int),
+            ("ins", int),
+            ("del", int),
+            ("consensus", "U1"),
+        ],
+    )
+
+    # Extract nucleotide counts for consensus calculation
+    nucleotides = np.vstack([data[base] for base in "ACGT"]).T
+    total_coverage = nucleotides.sum(axis=1)
+    max_counts = nucleotides.max(axis=1)
+
+    # Update consensus column where conditions are met
+    mask = np.divide(max_counts, total_coverage, where=total_coverage > 0) >= 0.7
+    data["consensus"][mask] = np.array(["A", "C", "G", "T"])[nucleotides[mask].argmax(axis=1)]
+
+    return data
+
+
+def write_csv(matrix: NDArray, prefix: str) -> None:
+    """
+    Write the matrix to a csv file
+
+    Args:
+        matrix: NumPy array containing the mpileup results
+        output: Path to the output file
+    """
+    header = ["Position", "Reference", "A", "C", "G", "T", "Insertions", "Deletions", "Consensus"]
+    with open(f"{prefix}.tsv", "w", newline="", encoding="utf-8") as file:
+        writer = csv.writer(file, delimiter="\t")
+        writer.writerow(header)
+        writer.writerows(matrix)
+
+
+def main():
+    args = parse_args()
+    logger.info("Starting mpileup processing")
+    matrix = process_mpileup(args.alignment, args.reference)
+    write_csv(matrix, args.prefix)
+    logger.info("Mpileup processing completed")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/conf/modules.config b/conf/modules.config
@@ -1264,14 +1264,26 @@ process {
                 ext.prefix = { "${meta.id}_${meta.previous_step}" } // DON'T CHANGE
                 publishDir = [
                     path: { "${meta.iteration}" == "variant-calling" ?
-                    "${params.outdir}/variants/mapping-info/picard/${meta.sample}" :
+                    "${params.outdir}/variants/mapping-info/metrics/picard/${meta.sample}" :
                     "${params.outdir}/assembly/polishing/iterations/${meta.step}/metrics/picard/" },
                     mode: params.publish_dir_mode,
                     pattern: '*{metrics,pdf}',
                     saveAs: { filename -> params.prefix || params.global_prefix  ? "${params.global_prefix}-$filename" : filename }
                 ]
             }
 
+            withName: CUSTOM_MPILEUP {
+                ext.prefix = { "${meta.id}_${meta.previous_step}.vcf" } // DON'T CHANGE
+                publishDir = [
+                    path: { "${meta.iteration}" == "variant-calling" ?
+                    "${params.outdir}/variants/mapping-info/custom-vcf/${meta.sample}" :
+                    "${params.outdir}/assembly/polishing/iterations/${meta.step}/custom-vcf/" },
+                    mode: params.publish_dir_mode,
+                    pattern: '*tsv',
+                    saveAs: { filename -> params.prefix || params.global_prefix  ? "${params.global_prefix}-$filename" : filename }
+                ]
+            }
+
             withName: ".*BAM_STATS_METRICS:BAM_STATS_SAMTOOLS:SAMTOOLS_STATS" {
                 ext.prefix = { "${meta.id}_${meta.previous_step}" } // DON'T CHANGE
                 publishDir = [
@@ -1471,9 +1483,7 @@ process {
             ext.args2 = 10
             ext.prefix = { "${meta.id}_${meta.previous_step}" } // DON'T CHANGE
             publishDir = [
-                path: { "${meta.iteration}" == "variant-calling" ?
-                    "${params.outdir}/consensus/mask/${meta.sample}" :
-                    "${params.outdir}/assembly/polishing/iterations/${meta.step}/consensus/mask" },
+                path: { "${params.outdir}/consensus/mask/${meta.step}/${meta.sample}"},
                 mode: params.publish_dir_mode,
                 enabled: params.save_intermediate_polishing,
                 pattern: '*.{mpileup}',
@@ -1484,9 +1494,7 @@ process {
         withName: BEDTOOLS_MERGE {
             ext.prefix = { "${meta.id}_${meta.previous_step}.merged" } // DON'T CHANGE
             publishDir = [
-                path: { "${meta.iteration}" == "variant-calling" ?
-                    "${params.outdir}/consensus/mask/${meta.sample}" :
-                    "${params.outdir}/assembly/polishing/iterations/${meta.step}/consensus/mask" },
+                path: { "${params.outdir}/consensus/mask/${meta.step}/${meta.sample}"},
                 mode: params.publish_dir_mode,
                 enabled: params.save_intermediate_polishing,
                 pattern: '*.{bed}',
@@ -1497,9 +1505,7 @@ process {
         withName: BEDTOOLS_MASKFASTA {
             ext.prefix = { "${meta.id}_${meta.previous_step}.masked" } // DON'T CHANGE
             publishDir = [
-                path: { "${meta.iteration}" == "variant-calling" ?
-                    "${params.outdir}/consensus/mask/${meta.sample}" :
-                    "${params.outdir}/assembly/polishing/iterations/${meta.step}/consensus/mask" },
+                path: { "${params.outdir}/consensus/mask/${meta.step}/${meta.sample}"},
                 mode: params.publish_dir_mode,
                 enabled: false,
                 pattern: '*.{fa}',
@@ -1538,24 +1544,11 @@ process {
             ].join(' ').trim()
             ext.prefix = { "${meta.id}_it${meta.iteration}" } // DON'T CHANGE
             publishDir = [
-                [
-                    path: { "${meta.iteration}" == "variant-calling" ?
-                    "${params.outdir}/consensus/mask/${meta.sample}" :
-                    "${params.outdir}/assembly/polishing/iterations/${meta.step}/consensus/mask/${meta.sample}" },
-                    mode: params.publish_dir_mode,
-                    enabled: params.save_intermediate_polishing,
-                    pattern: "*.txt",
-                    saveAs: { filename -> params.prefix || params.global_prefix  ? "${params.global_prefix}-$filename" : filename }
-                ],
-                [
-                    path: { "${meta.iteration}" == "variant-calling" ?
-                        "${params.outdir}/consensus/mask/${meta.sample}" :
-                        "${params.outdir}/assembly/polishing/iterations/${meta.step}/consensus/mask/${meta.sample}" },
-                    mode: params.publish_dir_mode,
-                    enabled: params.save_intermediate_polishing,
-                    pattern: "*.mpileup",
-                    saveAs: { filename -> params.prefix || params.global_prefix  ? "${params.global_prefix}-$filename" : filename }
-                ]
+                path: { "${params.outdir}/consensus/mask/${meta.step}/${meta.sample}"},
+                mode: params.publish_dir_mode,
+                enabled: params.save_intermediate_polishing,
+                pattern: "*.{txt,mpileup}",
+                saveAs: { filename -> params.prefix || params.global_prefix  ? "${params.global_prefix}-$filename" : filename }
             ]
         }
 
@@ -1596,7 +1589,7 @@ process {
         }
 
         withName: QUAST_QC {
-            ext.args      = "--min-contig 0"
+            ext.args   = "--min-contig 0"
             publishDir = [
                 path: { "${params.outdir}/consensus/quality_control/quast/${meta.sample}/${meta.step}/" },
                 mode: params.publish_dir_mode,
diff --git a/docs/output.md b/docs/output.md
@@ -574,6 +574,21 @@ __Summary statistics__:
         - `*.CollectMultipleMetrics.*`: Alignment QC files from picard CollectMultipleMetrics in `*_metrics` textual format.
         - `*.pdf` plots for metrics obtained from CollectMultipleMetrics.
 
+#### Custom - mpileup like file
+
+To facilitate the intra host analysis, a mpileup like file is generated. This file contains the depth of every nucletoride at each position of the reference.
+
+???- abstract "Output files - variants"
+
+    - `variants/mapping-info/custom-vcf/<sample-id>`
+        - `*.tsv`: A custom tsv file containing the depth of every nucleotide at each position of the reference.
+
+
+???- abstract "Output files - iterations"
+
+    - `assembly/polishing/iterations/it#/custom-vcf/<sample-id>`
+        - `*.tsv`: A custom tsv file containing the depth of every nucleotide at each position of the reference.
+
 #### Mosdepth - Coverage
 
 [mosdepth](https://github.com/brentp/mosdepth) is a fast BAM/CRAM depth calculation for WGS, exome, or targeted sequencing. mosdepth is used in this pipeline to obtain genome-wide coverage values in 200bp windows. The results are rendered in MultiQC (genome-wide coverage).
@@ -645,25 +660,15 @@ The consensus sequences are generated by [`BCFTools`](http://samtools.github.io/
 
 `BCFtools` will use the filtered variants file whereas, `iVar` will redetermine the variants to collapse in the consensus using their own workflow, read more about their differences in the [consensus calling section](./workflow/variant_and_refinement.md#consensus-calling).
 
-???- abstract "Output files - variants"
+???- abstract "Output files - iterations & variants"
 
     - `consensus`
-        - `seq/`
-            - `<sample-id>/<sample-id>_<constrain-id>-CONSTRAIN_itvariant_calling.fa`: A fasta file containing the consensus sequence.
-        - `mask/`
-            - `<sample-id>/<sample-id>_<constrain-id>-CONSTRAIN*.qual.txt`: A log file of the consensus run containing statistics. [`iVar` only]
-            - `<sample-id>/<sample-id>_<constrain-id>-CONSTRAIN*.bed`: A bed file containing the masked regions. [`BCFtools` only]
-            - `<sample-id>/<sample-id>_<constrain-id>-CONSTRAIN*.mpileup`: A mpileup file containing information on the depth and the quality of each alinged base.
-
-???- abstract "Output files - iterations"
-
-        - `assembly/polishing/iterations/it#/consensus`
-            - `seq/`
-                - `<sample-id>/<sample-id>_cl#_itvariant_calling.fa`: A fasta file containing the consensus sequence.
-            - `mask/`
-                - `<sample-id>/<sample-id>_cl#_it*.qual.txt`: A log file of the consensus run containing statistics. [`iVar` only]
-                - `<sample-id>/<sample-id>_cl#_it*.bed`: A bed file containing the masked regions. [`BCFtools` only]
-                - `<sample-id>/<sample-id>_cl#_it*.mpileup`: A mpileup file containing information on the depth and the quality of each alinged base.
+        - `seq/<it# | consensus | singleton | constrain>/ `
+            - `<sample-id>/*.fasta`: A fasta file containing the consensus sequence.
+        - `mask/<it# | consensus | singleton | constrain>`
+            - `<sample-id>/*.qual.txt`: A log file of the consensus run containing statistics. [`iVar` only]
+            - `<sample-id>/*.bed`: A bed file containing the masked regions. [`BCFtools` only]
+            - `<sample-id>/*.mpileup`: A mpileup file containing information on the depth and the quality of each alinged base.
 
 
 ## Consensus Quality control
diff --git a/modules/local/custom_mpileup/environment.yml b/modules/local/custom_mpileup/environment.yml
@@ -0,0 +1,9 @@
+name: custom_mpileup
+channels:
+  - conda-forge
+  - bioconda
+dependencies:
+  - bioconda::argparse
+  - bioconda::pysam
+  - bioconda::pysamstats
+  - conda-forge::pandas
diff --git a/modules/local/custom_mpileup/main.nf b/modules/local/custom_mpileup/main.nf
@@ -0,0 +1,52 @@
+process CUSTOM_MPILEUP {
+    tag "$meta.id"
+    label 'process_single'
+
+    conda "${moduleDir}/environment.yml"
+    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
+        'https://depot.galaxyproject.org/singularity/pysamstats:1.1.2--py39he47c912_12':
+        'quay.io/biocontainers/pysamstats:1.1.2--py311h0152c62_12' }"
+
+    input:
+    tuple val(meta), path(bam), path(ref)
+
+    output:
+    tuple val(meta), path("*.tsv"), emit: tsv
+    path "versions.yml"           , emit: versions
+
+    when:
+    task.ext.when == null || task.ext.when
+
+    script:
+    def args = task.ext.args ?: ''
+    def prefix = task.ext.prefix ?: "${meta.id}"
+    """
+    custom_mpileup.py \\
+        $args \\
+        --alignment ${bam} \\
+        --reference ${ref} \\
+        --prefix ${prefix}
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        python: \$(python --version | sed 's/Python //g')
+        pandas: \$(pip show pandas | grep Version: | sed 's/Version: //g')
+        pysam: \$(pip show pysam | grep Version: | sed 's/Version: //g')
+        pysamstats: \$(pip show pysamstats | grep Version: | sed 's/Version: //g')
+    END_VERSIONS
+    """
+
+    stub:
+    def args = task.ext.args ?: ''
+    def prefix = task.ext.prefix ?: "${meta.id}"
+    """
+    touch ${prefix}.tsv
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        python: \$(python --version | sed 's/Python //g')
+        pandas: \$(pip show pandas | grep Version: | sed 's/Version: //g')
+        pysam: \$(pip show pysam | grep Version: | sed 's/Version: //g')
+        pysamstats: \$(pip show pysamstats | grep Version: | sed 's/Version: //g')
+    """
+}
diff --git a/subworkflows/local/bam_stats_metrics.nf b/subworkflows/local/bam_stats_metrics.nf
@@ -2,6 +2,7 @@ include { SAMTOOLS_INDEX                } from '../../modules/nf-core/samtools/i
 include { PICARD_COLLECTMULTIPLEMETRICS } from '../../modules/nf-core/picard/collectmultiplemetrics/main'
 include { MOSDEPTH                      } from '../../modules/nf-core/mosdepth/main'
 include { BAM_STATS_SAMTOOLS            } from '../nf-core/bam_stats_samtools/main'
+include { CUSTOM_MPILEUP                } from '../../modules/local/custom_mpileup/main'
 
 workflow BAM_STATS_METRICS {
 
@@ -27,6 +28,8 @@ workflow BAM_STATS_METRICS {
             bam_bai_bed: [meta, bam, bai, []]
         }
 
+    CUSTOM_MPILEUP (sort_bam_ref)
+
     PICARD_COLLECTMULTIPLEMETRICS ( input_metrics.bam_bai, input_metrics.ref, [[:], []] )
     ch_versions  = ch_versions.mix(PICARD_COLLECTMULTIPLEMETRICS.out.versions)
 
