added keep, refactored mity call, normalise, util to use classes

Author: trentzz

Dockerfile.dev

@@ -6,26 +6,26 @@ ARG TAG
 # Install freebayes dependencies
 RUN apt-get update -yqq && \
     apt-get install -yqq \
-        make build-essential libssl-dev zlib1g-dev libbz2-dev \
-        libreadline-dev libsqlite3-dev wget curl llvm libncurses5-dev libncursesw5-dev \
-        xz-utils tk-dev libffi-dev liblzma-dev python3-openssl git tabix \
-        freebayes && \
+    make build-essential libssl-dev zlib1g-dev libbz2-dev \
+    libreadline-dev libsqlite3-dev wget curl llvm libncurses5-dev libncursesw5-dev \
+    xz-utils tk-dev libffi-dev liblzma-dev python3-openssl git tabix \
+    freebayes && \
     apt-get clean
 
 # Install gsort
 RUN \
     cd /usr/local/bin && \
     curl -s https://api.github.com/repos/brentp/gsort/releases/latest \
-      | grep browser_download_url \
-      | grep -i $(uname) \
-      | cut -d '"' -f 4 \
-      | wget -O gsort -qi - && \
+    | grep browser_download_url \
+    | grep -i $(uname) \
+    | cut -d '"' -f 4 \
+    | wget -O gsort -qi - && \
     chmod +x gsort
 
 # Install mity from dev server (but first install the previous version to get the dependencies from pypi)
 # RUN pip install mitywgs==0.4.0
 # mitywgs 0.4.0 requires pyvcf which doesn't allow testing of python 3.10.6 so here I'm installing dependencies manually
-RUN pip install pandas xlsxwriter pyfastx scipy pysam
+RUN pip install pandas xlsxwriter pyfastx scipy pysam pyyaml
 RUN pip install -i https://test.pypi.org/simple/ mitywgs-test==${TAG}
 
 WORKDIR /home

mitylib/call.py

@@ -157,7 +157,7 @@ def set_mity_cmd(self):
 
         mity_cmd += " " + " ".join(self.files)
         mity_cmd += '"'
-        mity_cmd = mity_cmd.replace("/", "\/")
+        mity_cmd = mity_cmd.replace("/", "\\/")
 
         logger.debug(mity_cmd)
 
@@ -215,26 +215,32 @@ def run_checks(self):
 
     def bam_has_rg(self, bam):
         """
-        Does the BAM or CRAM File have an @RG header? This is critical for mity
-        to correctly call variants.
+        Check whether a BAM or CRAM file contains a valid @RG header, which is critical for accurate variant calling with mity.
 
-        :param bam: str: path to bam or cram file
-        :return: True/False
-        >>> bam_has_RG('NA12878.alt_bwamem_GRCh38DH.20150718.CEU.low_coverage.chrM.bam')
+        Parameters:
+            - bam (str): Path to a BAM or CRAM file.
+
+        Returns:
+            - bool: True if the file has a valid @RG header, False otherwise.
         """
         r = pysam.AlignmentFile(bam, "rb")
         return len(r.header["RG"]) > 0
 
     def bam_get_mt_contig(self, bam, as_string=False):
         """
-        get the mitochondrial contig name and length from a BAM file
-        :param bam: path to a bam or cram file
-        :return: a tuple of contig name as str and length as int
-
-        >>> bam_get_mt_contig('NA12878.alt_bwamem_GRCh38DH.20150718.CEU.low_coverage.chrM.bam', False)
-        ('chrM', 16569)
-        >>> bam_get_mt_contig('NA12878.alt_bwamem_GRCh38DH.20150718.CEU.low_coverage.chrM.bam', True)
-        'chrM:1-16569'
+        Retrieve mitochondrial contig information from a BAM or CRAM file.
+
+        Parameters:
+            - bam (str): Path to a BAM or CRAM file.
+
+        Keyword Arguments:
+            - with_coordinates (bool): If True, the result includes the contig
+            name with coordinates (e.g., 'chrM:1-16569'). If False, it returns a
+            tuple containing the contig name (str) and its length (int).
+
+        Returns:
+            - If with_coordinates is False, a tuple (contig_name, contig_length).
+            - If with_coordinates is True, a string with coordinates.
         """
         r = pysam.AlignmentFile(bam, "rb")
         chroms = [str(record.get("SN")) for record in r.header["SQ"]]

mitylib/commands.py

@@ -4,7 +4,7 @@
 
 
 Usage: See the online manual for details: http://github.com/KCCG/mity
-Authors: Clare Puttick, Mark Cowley
+Authors: Clare Puttick, Mark Cowley, Trent Zeng, Christian Fares
 License: MIT
 """
 import argparse
@@ -13,13 +13,6 @@
 from mitylib import call, normalise, report, merge, util
 from ._version import __version__
 
-__all__ = []
-
-
-def public(fn):
-    __all__.append(fn.__name__)
-    return fn
-
 
 usage = __doc__.split("\n\n\n", maxsplit=1)
 usage[-1] += "Version: " + __version__
@@ -34,8 +27,6 @@ def public(fn):
 
 # call -------------------------------------------------------------------------
 
-do_call = public(call.do_call)
-
 
 def _cmd_call(args):
     """Call mitochondrial variants"""
@@ -157,7 +148,8 @@ def _cmd_call(args):
     action="store",
     type=str,
     default=None,
-    help="A text file of BAM files to be processed. The path to each file should be on one row per Region of MT genome to call variants in. "
+    help="A text file of BAM files to be processed. The path to each file should be"
+    "on one row per Region of MT genome to call variants in. "
     "If unset will call variants in entire MT genome as specified in BAM header. "
     "Default: Entire MT genome. ",
     dest="region",
@@ -166,8 +158,6 @@ def _cmd_call(args):
 
 # normalise --------------------------------------------------------------------
 
-do_normalise = public(normalise.do_normalise)
-
 
 def _cmd_normalise(args):
     """Normalise & FILTER mitochondrial variants"""
@@ -177,13 +167,14 @@ def _cmd_normalise(args):
     genome = util.MityUtil.select_reference_genome(args.reference, None)
     args.reference = util.MityUtil.select_reference_fasta(args.reference, None)
 
-    normalise.do_normalise(
+    normalise.Normalise(
         args.debug,
         args.vcf,
         args.reference,
         genome,
         args.outfile,
         args.allsamples,
+        args.keep,
         p=args.p,
     )
 
@@ -206,6 +197,13 @@ def _cmd_normalise(args):
     required=False,
     help="PASS in the filter requires all samples to pass instead of just one",
 )
+P_normalise.add_argument(
+    "-k",
+    "--keep",
+    action="store_true",
+    required=False,
+    help="Keep all intermediate files",
+)
 P_normalise.add_argument(
     "--p",
     action="store",
@@ -227,15 +225,13 @@ def _cmd_normalise(args):
 
 # report -----------------------------------------------------------------------
 
-do_report = public(report.do_report)
-
 
 def _cmd_report(args):
     """Generate mity report"""
     logging.info("mity %s", __version__)
     logging.info("Generating mity report")
-    report.do_report(
-        args.debug, args.vcf, args.prefix, args.min_vaf, args.out_folder_path
+    report.Report(
+        args.debug, args.vcf, args.prefix, args.min_vaf, args.out_folder_path, args.keep
     )
 
 
@@ -265,13 +261,18 @@ def _cmd_report(args):
 P_report.add_argument(
     "vcf", action="append", nargs="+", help="mity vcf files to create a report from"
 )
+P_report.add_argument(
+    "-k",
+    "--keep",
+    action="store_true",
+    required=False,
+    help="Keep all intermediate files",
+)
 P_report.set_defaults(func=_cmd_report)
 
 
 # merge -----------------------------------------------------------------------
 
-do_merge = public(merge.do_merge)
-
 
 def _cmd_merge(args):
     """Merging mity VCF with nuclear VCF"""

mitylib/mity

@@ -9,15 +9,17 @@ License: MIT
 """
 
 import logging
-from mitylib import commands
-from mitylib import util
 import sys
+from mitylib import commands
 
 assert sys.version_info >= (3, 7, 4)
 
+
 def main():
+    """
+    Entry function for mity.
+    """
     logging.basicConfig(level=logging.INFO, format="%(message)s")
-    util.check_dependencies()
 
     args = commands.parse_args()
 
@@ -27,5 +29,5 @@ def main():
         args.func(args)
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()

mitylib/normalise.py

@@ -2,14 +2,13 @@
 
 import subprocess
 import logging
-from .util import MityUtil
-from scipy.stats import binom
+import os.path
 from math import isinf
+from scipy.stats import binom
 import numpy as np
-import configparser
 import pysam
 import pysam.bcftools
-import os.path
+from mitylib.util import MityUtil
 
 
 # LOGGING
@@ -37,6 +36,7 @@ def __init__(
         genome,
         output_file=None,
         allsamples=False,
+        keep=False,
         p=P_VAL,
     ):
         self.debug = debug
@@ -45,6 +45,7 @@ def __init__(
         self.genome = genome
         self.output_file = output_file
         self.allsamples = allsamples
+        self.keep = keep
         self.p = p
 
         self.normalised_vcf_obj = None
@@ -88,9 +89,10 @@ def run(self):
         self.do_gsort()
         MityUtil.tabix(output_file)
 
-        # remove temporary files
-        os.remove(self.filtered_vcf_name)
-        os.remove(self.normalised_vcf_name)
+        if not self.keep:
+            # remove temporary files
+            os.remove(self.filtered_vcf_name)
+            os.remove(self.normalised_vcf_name)
 
     def run_bcftools_norm(self):
         """
@@ -210,25 +212,14 @@ def add_headers(self):
 
     def mity_qual(self, AO, DP, p):
         """
-        Compute variant quality
-        :param AO: (int) number of alternative reads
-        :param DP: (int) total read depth
-        :return: (float) phred-scaled quality score
-
-        >>> [mity_qual(x,10) for x in [1,2,3,4,5,6,7,8,9,10]]
-        '[37.49, 60.22, 84.78, 110.97, 138.75, 168.16, 199.4, 232.92, 269.9, 296.89]'
-        >>> [mity_qual(x,100) for x in [5,10,15,20,25,30,35,40,45,50]]
-        '[71.87, 156.08, 251.23, 354.34, 463.9, 579.05, 699.21, 824.04, 953.32, 1086.94]'
-        >>> [mity_qual(x,1000) for x in [5,10,15,20,25,30,35,40,45,50]]
-        '[17.84, 50.97, 93.59, 142.89, 197.36, 256.02, 318.25, 383.57, 451.61, 522.1]'
-        >>> [mity_qual(x,10000) for x in [5,10,15,20,25,30,35,40,45,50]]
-        '[0.0, 0.05, 0.74, 3.56, 9.51, 18.73, 31.0, 46.04, 63.57, 83.37]'
-        >>> mity_qual(118,118)
-        '3211.76'
-        >>> mity_qual(119,119)
-        '3220'
-        >>> mity_qual(10000,10000)
-        '3220'
+        Calculate the Phred-scaled quality score for a genetic variant.
+
+        Parameters:
+            - AO (int): Number of alternative reads.
+            - DP (int): Total read depth.
+
+        Returns:
+            - float: Phred-scaled quality score.
         """
         q = 0.0
         AO = int(AO)
@@ -248,6 +239,15 @@ def mity_qual(self, AO, DP, p):
         return q
 
     def add_filter(self, variant):
+        """
+        Adds filter to variant, and filters for samples.
+
+        Parameters:
+            - variant (VariantRecord): Variant to add filters to
+
+        Returns:
+            - VariantRecord: Variant with filters added
+        """
         # filter dictionary
         # each value represents the number of samples that pass each field
         # starts with all samples passing
@@ -391,5 +391,5 @@ def do_gsort(self):
         """
         Run gsort.
         """
-        gsort_cmd = f"gsort {self.filtered_vcf_name} {self.genome} | bgzip -cf > {self.gsorted_name}"
+        gsort_cmd = f"gsort {self.filtered_vcf_name} {self.genome} | bgzip -cf > {self.output_file}"
         subprocess.run(gsort_cmd, shell=True, check=False)