Merge pull request #52 from KCCG/release/1.1.0

Release/1.1.0

Author: christianfares

Dockerfile

@@ -3,12 +3,12 @@ FROM python:3.10.6-slim
 
 ARG TAG
 
-# Install freebayes dependencies
+# Install 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 \
+        xz-utils tk-dev libffi-dev liblzma-dev python3-openssl git tabix vcfanno \
         freebayes && \
     apt-get clean
 

Dockerfile.dev

@@ -3,12 +3,12 @@ FROM python:3.10.6-slim
 
 ARG TAG
 
-# Install freebayes dependencies
+# Install 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 \
+    xz-utils tk-dev libffi-dev liblzma-dev python3-openssl git tabix vcfanno \
     freebayes && \
     apt-get clean
 

mitylib/annot_chrm/chrm_anticodon_positions.bed.gz

@@ -1,4 +1,5 @@
 """Mitochondrial variant calling."""
+
 import subprocess
 import logging
 import os.path
@@ -78,7 +79,7 @@ def run(self):
         else:
             logger.setLevel(logging.INFO)
 
-        if (self.bam_list):
+        if self.bam_list:
             self.get_files_from_list()
         self.run_checks()
         self.set_strings()
@@ -228,7 +229,8 @@ def run_checks(self):
 
     def bam_has_rg(self, bam):
         """
-        Check whether a BAM or CRAM file contains a valid @RG header, which is critical for accurate variant calling with mity.
+        Check whether a BAM or CRAM file contains a valid @RG header,
+        which is critical for accurate variant calling with mity.
 
         Parameters:
             - bam (str): Path to a BAM or CRAM file.
@@ -275,7 +277,8 @@ def make_prefix(self, file_name: str, prefix: str = None):
 
         Parameters:
             file_name (str): The filename, including extensions (e.g., .vcf, .bam, .cram, .bed).
-            prefix (str, optional): An optional custom prefix. If None, the function generates a prefix from the file name.
+            prefix (str, optional): An optional custom prefix. If None, the function generates a
+                                    prefix from the file name.
 
         Returns:
             str: The generated or custom prefix for the Mity function.
@@ -291,24 +294,25 @@ def make_prefix(self, file_name: str, prefix: str = None):
             )
         else:
             raise ValueError("Unsupported file type")
-        
+
     def get_files_from_list(self, die: bool = True):
         """
         Get the list of BAM / CRAM files from the provided file
         """
         if len(self.files) > 1:
-            raise ValueError("--bam-file-list Argument expects only 1 file to be provided.")
-        with open(self.files[0], 'r') as f:
+            raise ValueError(
+                "--bam-file-list Argument expects only 1 file to be provided."
+            )
+        with open(self.files[0], "r") as f:
             self.files = f.read().splitlines()
 
-
     def check_missing_file(self, file_list: str, die: bool = True):
         """
         Check if input files exist.
         """
         missing_files = []
         for item in file_list:
-            if item.lower().startswith('http'):
+            if item.lower().startswith("http"):
                 try:
                     urllib.request.urlopen(item).getcode()
                 except:

mitylib/annot_chrm/chrm_anticodon_positions.bed.gz.tbi

@@ -7,8 +7,10 @@
 Authors: Clare Puttick, Mark Cowley, Trent Zeng, Christian Fares
 License: MIT
 """
+
 import argparse
 import logging
+import os
 
 from mitylib import call, normalise, report, merge, util
 from ._version import __version__
@@ -247,10 +249,13 @@ def _cmd_report(args):
     report.Report(
         debug=args.debug,
         vcfs=args.vcf,
+        contig=args.contig,
         prefix=args.prefix,
         min_vaf=args.min_vaf,
         output_dir=args.output_dir,
         keep=args.keep,
+        vcfanno_config=args.vcfanno_config,
+        report_config=args.report_config,
     )
 
 
@@ -287,6 +292,25 @@ def _cmd_report(args):
     required=False,
     help="Keep all intermediate files",
 )
+P_report.add_argument(
+    "--contig",
+    choices=["MT", "chrM"],
+    default="MT",
+    required=False,
+    help="Contig used for annotation purposes",
+)
+P_report.add_argument(
+    "--custom-vcfanno-config",
+    action="store",
+    help="Provide a custom vcfanno-config.toml for custom annotations.",
+    dest="vcfanno_config",
+)
+P_report.add_argument(
+    "--custom-report-config",
+    action="store",
+    help="Provide a custom report-config.yaml for custom report generation.",
+    dest="report_config",
+)
 P_report.set_defaults(func=_cmd_report)
 
 
@@ -349,6 +373,201 @@ def _cmd_merge(args):
 )
 P_merge.set_defaults(func=_cmd_merge)
 
+
+# runall -----------------------------------------------------------------------
+
+
+def _cmd_runall(args):
+    """Run MITY call, normalise and report all in one go."""
+    logging.info("mity %s", __version__)
+    logging.info("mity runall")
+
+    genome = util.MityUtil.select_reference_genome(args.reference, None)
+    args.reference = util.MityUtil.select_reference_fasta(args.reference, None)
+
+    call.Call(
+        debug=args.debug,
+        files=args.files,
+        reference=args.reference,
+        genome=genome,
+        prefix=args.prefix,
+        min_mq=args.min_mq,
+        min_bq=args.min_bq,
+        min_af=args.min_af,
+        min_ac=args.min_ac,
+        p=args.p,
+        # normalise flag is set to true instead of running normalise separately
+        normalise=True,
+        output_dir=args.output_dir,
+        region=args.region,
+        bam_list=args.bam_file_list,
+        keep=args.keep,
+    )
+
+    logging.debug("mity call and normalise completed")
+
+    # This makes use of the uniform naming scheme of mity command outputs.
+    normalised_vcf_path = os.path.join(
+        args.output_dir, args.prefix + ".normalise.vcf.gz"
+    )
+
+    logging.debug("assumed mity normalise vcf output path is: %s", normalised_vcf_path)
+
+    # matching argparse quirk
+    normalised_vcf_path = [normalised_vcf_path]
+
+    report.Report(
+        debug=args.debug,
+        vcfs=normalised_vcf_path,
+        contig=args.contig,
+        prefix=args.prefix,
+        min_vaf=args.min_vaf,
+        output_dir=args.output_dir,
+        keep=args.keep,
+        vcfanno_config=args.vcfanno_config,
+        report_config=args.report_config,
+    )
+
+
+P_runall = AP_subparsers.add_parser("runall", help=_cmd_runall.__doc__)
+P_runall.add_argument(
+    "-d", "--debug", action="store_true", help="Enter debug mode", required=False
+)
+P_runall.add_argument(
+    "files",
+    action="append",
+    nargs="+",
+    help="BAM / CRAM files to run the analysis on. If --bam-file-list is included, this argument is the file containing the list of bam/cram files.",
+)
+P_runall.add_argument(
+    "--reference",
+    choices=["hs37d5", "hg19", "hg38", "mm10"],
+    default="hs37d5",
+    required=False,
+    help="Reference genome version to use. Default: hs37d5",
+)
+# For the runall command, we mandate that the prefix option is set. This is not
+# true for regular mity call, normalise or report separately.
+P_runall.add_argument(
+    "--prefix",
+    action="store",
+    required=True,
+    help="Output files will be named with PREFIX",
+)
+P_runall.add_argument(
+    "--min-mapping-quality",
+    action="store",
+    type=int,
+    default=30,
+    help="Exclude alignments from analysis if they have a "
+    "mapping quality less than MIN_MAPPING_QUALITY. "
+    "Default: 30",
+    dest="min_mq",
+)
+P_runall.add_argument(
+    "--min-base-quality",
+    action="store",
+    type=int,
+    default=24,
+    help="Exclude alleles from analysis if their supporting "
+    "base quality is less than MIN_BASE_QUALITY. "
+    "Default: 24",
+    dest="min_bq",
+)
+P_runall.add_argument(
+    "--min-alternate-fraction",
+    action="store",
+    type=float,
+    default=0.01,
+    help="Require at least MIN_ALTERNATE_FRACTION "
+    "observations supporting an alternate allele within "
+    "a single individual in the in order to evaluate the "
+    "position. Default: 0.01, range = [0,1]",
+    dest="min_af",
+)
+P_runall.add_argument(
+    "--min-alternate-count",
+    action="store",
+    type=int,
+    default=4,
+    help="Require at least MIN_ALTERNATE_COUNT observations "
+    "supporting an alternate allele within a single "
+    "individual in order to evaluate the position. "
+    "Default: 4",
+    dest="min_ac",
+)
+P_runall.add_argument(
+    "--p",
+    action="store",
+    type=float,
+    default=0.002,
+    help="Minimum noise level. This is used to calculate QUAL score. "
+    "Default: 0.002, range = [0,1]",
+    dest="p",
+)
+P_runall.add_argument(
+    "--output-dir",
+    action="store",
+    type=str,
+    default=".",
+    help="Output files will be saved in OUTPUT_DIR. " "Default: '.' ",
+    dest="output_dir",
+)
+P_runall.add_argument(
+    "--region",
+    action="store",
+    type=str,
+    default=None,
+    help="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",
+)
+P_runall.add_argument(
+    "--bam-file-list",
+    action="store_true",
+    default=False,
+    help="Treat the file as a text file of BAM files to be processed."
+    " The path to each file should be on one row per bam file.",
+    dest="bam_file_list",
+)
+P_runall.add_argument(
+    "-k",
+    "--keep",
+    action="store_true",
+    required=False,
+    help="Keep all intermediate files",
+)
+P_runall.add_argument(
+    "--min_vaf",
+    action="store",
+    type=float,
+    default=0,
+    help="A variant must have at least this VAF to be included in the report. Default: "
+    "0.",
+)
+P_runall.add_argument(
+    "--contig",
+    choices=["MT", "chrM"],
+    default="MT",
+    required=False,
+    help="Contig used for annotation purposes",
+)
+P_runall.add_argument(
+    "--custom-vcfanno-config",
+    action="store",
+    help="Provide a custom vcfanno-config.toml for custom annotations.",
+    dest="vcfanno_config",
+)
+P_runall.add_argument(
+    "--custom-report-config",
+    action="store",
+    help="Provide a custom report-config.yaml for custom report generation.",
+    dest="report_config",
+)
+P_runall.set_defaults(func=_cmd_runall)
+
+
 # version ----------------------------------------------------------------------
 
 

mitylib/annot_chrm/chrm_gtf_annotations.bed.gz

@@ -0,0 +1,37 @@
+[[annotation]]
+file="annot_chrm/chrm_anticodon_positions.bed.gz"
+columns=[4]
+names=["anticodon"]
+ops=["self"]
+
+[[annotation]]
+file="annot_chrm/chrm_gtf_annotations.bed.gz"
+columns=[4,5]
+names=["GENE","GENE_BIOTYPE"]
+ops=["self", "self"]
+
+[[annotation]]
+file="annot_chrm/chrm_haplotype_data.vcf.gz"
+fields=["phylotree_mut","phylotree_haplotype"]
+ops=["self", "self"]
+
+[[annotation]]
+file="annot_chrm/chrm_mgrb_variants.vcf.gz"
+fields = ["MGRB_FILTER","MGRB_AN","MGRB_AC","MGRB_frequency"]
+ops = ["self","self","self","self"]
+
+[[annotation]]
+file="annot_chrm/chrm_mito_dna_func_loc.bed.gz"
+columns=[4,7]
+names=["Map_Locus", "Description"]
+ops=["self", "self"]
+
+[[annotation]]
+file="annot_chrm/chrm_mitomap_panel_annotations.vcf.gz"
+fields = ["locus_mitomap","num_references_mitomap","variant_amino_acid_change_mitomap","codon_position_mitomap","codon_number_mitomap","disease_mitomap","num_disease_references_mitomap","RNA_mitomap","homoplasmy_mitomap","heteroplasmy_mitomap","status_mitomap","disease_amino_acid_change_mitomap","GenBank_frequency_mitomap","commercial_panels"]
+ops=["self","self","self","self","self","self","self","self","self","self","self","self","self","self"]
+
+[[annotation]]
+file="annot_chrm/chrm_mitotip_score_fixed_del.vcf.gz"
+fields=["MitoTip_score","MitoTip_percentile","MitoTip_interpretation"]
+ops=["self","self","self"]