Develop

.gitignore

@@ -1,6 +1,31 @@
 # Ignore filetypes
 *.pyc
+*.pyo
+*.pyd
+__pycache__/
+
+# Virtual environments
 /python2env/
+/.venv/
+/venv/
+/env/
+
+# IDEs
 /.ipynb_checkpoints/
 /.vscode/
-/.idea/
+/.idea/
+
+# Test & coverage artifacts
+.coverage
+.coverage.*
+htmlcov/
+.pytest_cache/
+
+# NEAT log files
+*.log
+
+# Build / packaging
+dist/
+build/
+*.egg-info/
+*.egg

README.md

@@ -1,12 +1,12 @@
-# The NEAT Project v4.3.6
+# The NEAT Project v4.4
 
-Welcome to the NEAT project, the NExt-generation sequencing Analysis Toolkit, version 4.3.6. This release of NEAT 4.3.5 includes several fixes and a little bit of restructuring, including a parallel process for running `neat read-simulator`. Our tests show much improved performance. If the logs seem excessive, you might try using the `--log-level ERROR` to reduce the output from the logs. See the [ChangeLog](ChangeLog.md) for notes. NEAT 4.3.5 is the official release of NEAT 4.0. It represents a lot of hard work from several contributors at NCSA and beyond. With the addition of parallel processing, we feel that the code is ready for production, and future releases will focus on compatibility, bug fixes, and testing. Future releases for the time being will be enumerations of 4.3.X.
+Welcome to the NEAT project, the NExt-generation sequencing Analysis Toolkit, version 4.4. NEAT 4.4 is the official release of NEAT 4.0. It represents a lot of hard work from several contributors at NCSA and beyond. With the addition of parallel processing, we feel that the code is ready for production, and future releases will focus on compatibility, bug fixes, and testing. Future releases for the time being will be enumerations of 4.4.X.
 
-## NEAT v4.3.5
+## NEAT v4.4
 
-NEAT 4.3.5 marked the officially 'complete' version of NEAT 4.3, implementing parallelization. To add parallelization to your run, simply add the `threads` parameter in your configuration file and run `read-simulator` as normal. NEAT will take care of the rest. You can customize the parameters in your configuration file, as needed.
+NEAT 4.4 fixes a few bugs related to NEAT. This release of NEAT includes several fixes to bugs and a little bit of restructuring, including more options for quality score modeling, a bacterial genome wrapper, additional tests, and more. Our tests show much improved performance. If the logs seem excessive, you might try using the `--log-level ERROR` to reduce the output from the logs. See the [ChangeLog](ChangeLog.md) for notes. 
 
-We have completed major revisions on NEAT since 3.4 and consider NEAT 4.3.5 to be a stable release, in that we will continue to update and provide bug fixes and support. We will consider new features and pull requests. Please include justification for major changes. See [contribute](CONTRIBUTING.md) for more information. If you'd like to use some of our code in your own, no problem! Just review the [license](LICENSE.md), first.
+We have completed major revisions on NEAT since 3.4 and consider NEAT 4.4 to be a stable release, in that we will continue to update and provide bug fixes and support. We will consider new features and pull requests. Please include justification for major changes. See [contribute](CONTRIBUTING.md) for more information. If you'd like to use some of our code in your own, no problem! Just review the [license](LICENSE.md), first.
 
 We've deprecated NEAT's command-line interface options for the most part, opting to simplify things with configuration files. If you require the CLI for legacy purposes, NEAT 3.4 was our last release to be fully supported via command-line interface. Please convert your CLI commands to the corresponding configuration file for future runs.
 
@@ -22,8 +22,8 @@ To cite this work, please use:
 
 ## Table of Contents
 
-* [The NEAT Project v4.3.6](#the-neat-project-v436)
-* [NEAT v4.3.5](#neat-v435)
+* [The NEAT Project v4.4](#the-neat-project-v436)
+* [NEAT v4.4](#neat-v435)
 * [Table of Contents](#table-of-contents)
   * [Prerequisites](#prerequisites)
   * [Installation](#installation)
@@ -246,7 +246,7 @@ Features:
 
 ### Estimated runtimes
 
-To give users a sense of how long `neat read-simulator` runs may take, we benchmarked NEAT 4.3.5 on several reference genomes. All runs were paired-end, with read length of 150 bp, coverage of 10, fragment mean of 300 bp, and fragment standard deviation of 50 bp. Runtimes are reported as the average across three unique runs (`Avg. time (ms)`) and the corresponding runtime in minutes. Cells marked with N/A indicate that NEAT was not able to run to completion.
+To give users a sense of how long `neat read-simulator` runs may take, we benchmarked NEAT 4.4 on several reference genomes. All runs were paired-end, with read length of 150 bp, coverage of 10, fragment mean of 300 bp, and fragment standard deviation of 50 bp. Runtimes are reported as the average across three unique runs (`Avg. time (ms)`) and the corresponding runtime in minutes. Cells marked with N/A indicate that NEAT was not able to run to completion.
 
 Benchmarks were run on a System76 Meerkat with a 13th Gen Intel Core i3-1315U (8 logical cores, up to 4.50 GHz) and 16 GiB RAM, using a 512 GB SSD and Ubuntu 24.04.3 LTS (Linux kernel 6.14). Actual runtimes will vary depending on your hardware.
 
@@ -440,7 +440,7 @@ neat gen-mut-model reference.fa input_variants.vcf   \
         -o /home/me/models
 ```
 
-Trinucleotides are identified in the reference genome and the variant file. The mutation model uses trinucleotide context and selects mutation sites and alternate alleles with a transition matrix. Frequencies of each trinucleotide transition are calculated and output as a pickle file. Mutations are simulated to reflect the same context-dependent biases as the training data. In NEAT 4.3.5, we have only made minor optimizations to improve the speed, and the underlying statistical models are similar to those described in the original NEAT manuscript.
+Trinucleotides are identified in the reference genome and the variant file. The mutation model uses trinucleotide context and selects mutation sites and alternate alleles with a transition matrix. Frequencies of each trinucleotide transition are calculated and output as a pickle file. Mutations are simulated to reflect the same context-dependent biases as the training data. As of NEAT 4.4, we have only made minor optimizations to improve the speed, and the underlying statistical models are similar to those described in the original NEAT manuscript.
 
 | Option          | Description                                                                   |
 |-----------------|-------------------------------------------------------------------------------|
@@ -510,7 +510,7 @@ Finally, `-o` is the output directory for the model file and `-p` is the prefix
 
 ### `neat vcf_compare`
 
-Tool for comparing VCF files (Not yet implemented in NEAT 4.3.5).
+Tool for comparing VCF files (Not yet implemented in NEAT 4.4).
 
 ```bash
 neat vcf_compare
@@ -546,4 +546,4 @@ pytest -q tests
 Please see `CONTRIBUTING.md` for more information and further instructions.
 
 ### Note on Sensitive Patient Data
-ICGC's "Access Controlled Data" documentation can be found at <a href = https://docs.icgc.org/portal/access/ target="_blank">https://docs.icgc.org/portal/access/</a>. To have access to controlled germline data, a DACO must be submitted. Open tier data can be obtained without a DACO, but germline alleles that do not match the reference genome are masked and replaced with the reference allele. Controlled data includes unmasked germline alleles.
+ICGC's "Access Controlled Data" documentation can be found at <a href = https://docs.icgc.org/portal/access/ target="_blank">https://docs.icgc.org/portal/access/</a>. To have access to controlled germline data, a DACO must be submitted. Open tier data can be obtained without a DACO, but germline alleles that do not match the reference genome are masked and replaced with the reference allele. Controlled data includes unmasked germline alleles.

config_template/template_neat_config.yml

@@ -133,8 +133,7 @@ rng_seed: .
 # type: int | required = no
 min_mutations: .
 
-# Overwrite the output files, if they are named the same as the current run.
-# Default is to quit if files already exist to avoid data destruction
+# Overwrite output files if they already exist
 # type: bool | required = no | default = false
 overwrite_output: .
 
@@ -153,12 +152,12 @@ parallel_block_size: .
 threads: .
 
 # Delete the 'splits' directory after stitching completes
-# Note if threads == 1, this option has no effect.
+# Note: If threads == 1, this option has no effect.
 # type = bool | required: no | default = true
 cleanup_splits: .
 
 # Reuse existing files in '<out_dir>/splits' and skip the split step.
 # The directory must contain NEAT-generated files and must be in the output directory within "splits"
-# Note if threads == 1, this option has no effect.
+# Note: If threads == 1, this option has no effect.
 # type = bool | required: no | default = False
-reuse_splits: .
+reuse_splits: .

neat/bacterial_wrapper/__init__.py

@@ -0,0 +1,4 @@
+"""
+Load modules needed for other parts of the program
+"""
+from .runner import *

neat/bacterial_wrapper/runner.py

@@ -0,0 +1,203 @@
+import subprocess
+import gzip
+import shutil
+import yaml
+import pysam
+import unittest
+import os
+
+from pathlib import Path
+from typing import List
+from Bio import bgzf
+from Bio.bgzf import BgzfWriter, BgzfReader
+
+
+# Rearranges the bacterial chromosome by wrapping it around
+
+def wrapper(seq):
+    length = len(seq)
+
+    if (length % 2 == 0):
+        half_index = length // 2
+    else:
+        half_index = (length // 2) + 1
+
+    first_half = seq[:half_index]
+    second_half = seq[half_index:]
+
+    new_seq = second_half + first_half
+
+    return new_seq
+
+
+# Writes the newly rearranged chromosome's sequence to a new fasta file
+
+def write_fasta_file(new_seq, bacteria_name, fasta_header, output_dir_path, type):
+    fasta_file_name = f"{type}_{bacteria_name}.fna"
+    fasta_file_path = output_dir_path / fasta_file_name
+    fasta_file = open(fasta_file_path, "w")
+
+    fasta_file.write(fasta_header + "\n" + new_seq)
+
+    fasta_file.close()
+
+    return fasta_file_path
+
+
+# Writes a yml configuration file for the newly rearranged chromosome's fasta sequence
+# Splits the coverage in half for the reference and new config files
+# These use default values for all other parameters for NEAT
+
+def write_config_file(ref_config_file, rearranged_seq_file, orig_seq_file, bacteria_name, output_dir_path):
+    new_config_file_name = f"new_{bacteria_name}_config_test.yml"
+    old_config_file_name = f"{bacteria_name}_config_test.yml"
+
+    new_config_file_path = output_dir_path / new_config_file_name
+    old_config_file_path = output_dir_path / old_config_file_name
+
+    with open(ref_config_file, 'r') as ref_file, open(new_config_file_path, 'w') as new_file, open(old_config_file_path, 'w') as old_file:
+        for line in ref_file:
+            if line.find("reference:") != -1:
+                new_file.write(f"reference: {rearranged_seq_file}\n")
+                old_file.write(f"reference: {orig_seq_file}\n")
+            # elif line.find("coverage:") != -1:
+            #     if line.strip() == "coverage: .":
+            #         new_coverage = 5.0
+            #     else:
+            #         new_coverage = float((line.split(" "))[1].strip()) // 2
+
+            #     new_file.write(f"coverage: {new_coverage}\n")
+            #     old_file.write(f"coverage: {new_coverage}\n")
+            else:
+                new_file.write(line)
+                old_file.write(line)
+
+
+    ref_file.close()
+    new_file.close()
+    old_file.close()
+
+    return old_config_file_path, new_config_file_path
+
+
+# Runs the NEAT read simulator using the given config file
+
+def run_neat(config_file, output_dir, prefix):
+    subprocess.run(["neat", "read-simulator", "-c", config_file, "-o", output_dir + "/" + prefix])
+
+
+# General function for bacterial wrapper that calls all of the functions defined above
+
+def bacterial_wrapper(reference_file, bacteria_name, ref_config_file, output_dir):
+
+    orig_seq = ""
+
+    f = open(reference_file)
+    fasta_header = f.readline().strip()
+
+    plasmids = False
+
+    for line in f:
+        if line[0] != ">":
+            orig_seq += line.strip()
+        elif line.lower().find("plasmid") != -1:  # exclude plasmids from the sequence to be rearranged
+            plasmids = True
+            break
+
+    f.close()
+
+    output_dir_path = Path(output_dir)
+
+    rearranged_seq = wrapper(orig_seq)
+    rearranged_seq_file = write_fasta_file(rearranged_seq, bacteria_name, fasta_header, output_dir_path, "wrapped")
+
+    orig_seq_file = reference_file
+    if plasmids:
+        orig_seq_file = write_fasta_file(orig_seq, bacteria_name, fasta_header, output_dir_path, "orig")
+
+    config_files = write_config_file(ref_config_file, rearranged_seq_file, orig_seq_file, bacteria_name, output_dir_path)
+    old_config_file = config_files[0]
+    new_config_file = config_files[1]
+
+    run_neat(old_config_file, output_dir, "Regular")
+    run_neat(new_config_file, output_dir, "Wrapped")
+
+
+# Stitching all outputs together - Keshav's script
+
+def concat_fq(input_files: List[Path], dest: Path) -> None:
+
+    if not input_files:
+        # Nothing to do, and no error to throw
+        return
+
+    with gzip.open(dest, 'wt') as out_f:
+        for input_file in input_files:
+            with gzip.open(input_file, 'rt') as in_f:
+                shutil.copyfileobj(in_f, out_f)
+
+def merge_bam(bams: List[Path], dest: Path, threads: int) -> None:
+
+    if not bams:
+        return
+
+    unsorted = dest.with_suffix(".unsorted.bam")
+    pysam.merge("--no-PG", "-@", str(threads), "-f", str(unsorted), *map(str, bams))
+    pysam.sort("-@", str(threads), "-o", str(dest), str(unsorted))
+    unsorted.unlink(missing_ok=True)
+
+def merge_vcf(vcfs: List[Path], dest: Path) -> None:
+    if not vcfs:
+        return
+
+    first, *rest = vcfs
+    shutil.copy(first, dest)
+
+    with dest.open("ab") as out_f:
+        for vcf in rest:
+            with vcf.open("rb") as fh:
+                for line in fh:
+                    if not line.startswith(b"#"):
+                        out_f.write(line)
+
+def stitch_all_outputs(files: List[Path], output_dir) -> None:
+    fq1_list = []
+    fq2_list = []
+    vcf_list = []
+    bam_list = []
+
+    for file in files:
+        file_name = file.stem  # use stem to differentiate fq1 and fq2
+        suffixes = file.suffixes  # use suffixes to catch vcf and bam files
+
+        if "r2.fastq" in file_name:
+            fq2_list.append(file)
+        elif "r1.fastq" in file_name or ".fastq" in suffixes:
+            fq1_list.append(file)
+        elif ".vcf" in suffixes and ".tbi" not in suffixes:
+            vcf_list.append(file)
+        elif ".bam" in suffixes and ".bai" not in suffixes:
+            bam_list.append(file)
+
+    dest_fq1 = Path(f"{output_dir}/stitched_fq1.gz")
+    dest_bam = Path(f"{output_dir}/stitched.bam")
+    dest_vcf = Path(f"{output_dir}/stitched.vcf")
+
+    concat_fq(fq1_list, dest_fq1)
+
+    if (fq2_list):
+        dest_fq2 = Path(f"{output_dir}/stitched_fq2.gz")
+        concat_fq(fq2_list, dest_fq2)
+
+    merge_bam(bam_list, dest_bam, 2)
+    merge_vcf(vcf_list, dest_vcf)
+
+
+# Testing functions
+
+class TestWrapper(unittest.TestCase):
+    def test_even(self):
+        self.assertEqual(wrapper("ABBCBB"), "CBBABB")
+
+    def test_odd(self):
+        self.assertEqual(wrapper("ABBCBBC"), "BBCABBC")

neat/cli/commands/__init__.py

@@ -1,2 +1,3 @@
-"""Modules related to the subcommands' interfaces"""
+"""Modules related to the subcommands"""
+
 from .base import *