Modernize packaging, add test suite, remove shell-string system calls

.github/workflows/ci.yml

@@ -0,0 +1,33 @@
+name: CI
+
+on:
+  push:
+    branches: [master, main]
+  pull_request:
+
+jobs:
+  test:
+    name: Unit & smoke tests (Python ${{ matrix.python-version }})
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        python-version: ["3.9", "3.10", "3.11", "3.12"]
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Set up Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+
+      - name: Install package and test dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install -e .[test]
+
+      - name: Run unit and smoke tests
+        # Integration tests require samtools/bedtools/R/snakemake and are
+        # skipped automatically when those tools are not installed.
+        run: pytest -m "not integration" -v

MANIFEST.in

@@ -9,5 +9,7 @@ recursive-include bap *
 
 # Misc
 include Authors.rst
-include LICENSE.txt
+include LICENSE
 include NEWS
+include README.md
+include pyproject.toml

README.md

@@ -10,5 +10,27 @@ data files for downstream analyses for droplet-based single-cell ATAC-seq data.
 
 [See the bap wiki page](https://github.com/caleblareau/bap/wiki) for details and FAQs about bap.
 
+## Installation
+
+```
+pip install bap-atac          # from PyPI
+pip install -e .              # from a local checkout (development)
+```
+
+Requires Python >= 3.9. The pipeline also shells out to `samtools`, `bedtools`,
+`R`, and `snakemake`, which must be available on your `PATH`.
+
+## Testing
+
+Automated tests use `pytest`. Install the test extras and run the fast suite:
+
+```
+pip install -e .[test]
+pytest -m "not integration"   # unit + CLI smoke tests, no external tools needed
+pytest                        # also runs end-to-end tests when the external tools are present
+```
+
+See [`tests/README.md`](tests/README.md) for the full testing guide.
+
 
 

bap/bap2ProjectClass.py

@@ -8,8 +8,6 @@
 import itertools
 import time
 import platform
-from ruamel import yaml
-from pkg_resources import get_distribution
 from .bapHelp import *
 
 def getBfiles(bedtools_genome, blacklist_file, reference_genome, script_dir, supported_genomes):
@@ -75,7 +73,7 @@ def __init__(self, script_dir, supported_genomes, mode, input, output, name, nco
 		#----------------------------------
 		# Assign straightforward attributes
 		#----------------------------------
-		self.bap_version = get_distribution('bap-atac').version
+		self.bap_version = bap_version()
 		self.script_dir = script_dir
 		self.mode = mode
 		self.output = output
@@ -173,12 +171,12 @@ def __init__(self, script_dir, supported_genomes, mode, input, output, name, nco
 		else: 
 			click.echo(gettime() + "Could not identify this reference genome: %s" % self.reference_genome)
 			click.echo(gettime() + "Attempting to infer necessary input files from user specification.")
-			necessary = [bedtools_genome, blacklist_file, tss_file, macs2_genome_size, bs_genome]
+			necessary = [bedtools_genome, blacklist_file, tss_file]
 			if '' in necessary:
 				if reference_genome == '':
 					sys.exit("ERROR: specify valid reference genome with --reference-genome flag; QUITTING")
 				else:
-					sys.exit("ERROR: non-supported reference genome specified so these five must be validly specified: --bedtools-genome, --blacklist-file, --tss-file; QUITTING")
+					sys.exit("ERROR: non-supported reference genome specified so these three must be validly specified: --bedtools-genome, --blacklist-file, --tss-file; QUITTING")
 
 		if(reference_genome in ["hg19-mm10", "hg19_mm10_c"]):
 			self.speciesMix = "yes"

bap/bapFragProjectClass.py

@@ -8,8 +8,6 @@
 import itertools
 import time
 import platform
-from ruamel import yaml
-from pkg_resources import get_distribution
 from .bapHelp import *
 
 def getBfiles(bedtools_genome, blacklist_file, reference_genome, script_dir, supported_genomes):
@@ -69,7 +67,7 @@ def __init__(self, script_dir, supported_genomes, input, output, name, ncores, r
 		#----------------------------------
 		# Assign straightforward attributes
 		#----------------------------------
-		self.bap_version = get_distribution('bap-atac').version
+		self.bap_version = bap_version()
 		self.script_dir = script_dir
 		self.bamfile = input
 		self.name = name
@@ -202,7 +200,7 @@ def __init__(self, script_dir, supported_genomes, input, output, name, ncores, r
 		#----------------------------------
 		# Assign straightforward attributes
 		#----------------------------------
-		self.bap_version = get_distribution('bap').version
+		self.bap_version = bap_version()
 		self.script_dir = script_dir
 		self.bamfile = input
 		self.name = name

bap/bapHelp.py

@@ -1,23 +1,65 @@
-import itertools
 import time
 import shutil
 import re
 import os
 import sys
 import csv
 import gzip
-from ruamel import yaml
+import glob
+import subprocess
 from functools import partial
+from importlib.metadata import version, PackageNotFoundError
+
+
+def bap_version():
+	'''
+	Return the installed bap-atac version, or "unknown" if the package
+	metadata cannot be located (e.g. running from an uninstalled checkout).
+	'''
+	try:
+		return version("bap-atac")
+	except PackageNotFoundError:
+		return "unknown"
+
+
+def run_cmd(args, log_file=None, log_stderr_only=False, shell=False, check=True):
+	'''
+	Run an external command from a list of arguments (no shell by default),
+	raising a clear error if it fails. This replaces the historical
+	os.system(<string>) calls, avoiding shell-injection and quoting issues.
+
+	args:            list of command arguments (or a string when shell=True)
+	log_file:        optional path; command output is written here. By default
+	                 both stdout and stderr are captured (equivalent to `&>`).
+	log_stderr_only: when True, only stderr is captured to log_file and stdout
+	                 is inherited by the terminal (equivalent to `2>`).
+	shell:           only set True for commands that need a shell pipeline
+	check:           when True, exit with a clear message on a non-zero return
+	                 code; set False to let the caller inspect the result
+	'''
+	if log_file is not None:
+		with open(log_file, 'w') as fh:
+			if log_stderr_only:
+				completed = subprocess.run(args, stderr=fh, shell=shell)
+			else:
+				completed = subprocess.run(args, stdout=fh, stderr=subprocess.STDOUT, shell=shell)
+	else:
+		completed = subprocess.run(args, shell=shell)
+	if check and completed.returncode != 0:
+		printed = args if isinstance(args, str) else " ".join(str(a) for a in args)
+		sys.exit(gettime() + "ERROR: command failed (exit %d): %s" % (completed.returncode, printed))
+	return completed
+
 
 def string_hamming_distance(str1, str2):
 	'''
 	Fast hamming distance over 2 strings known to be of same length.
-	In information theory, the Hamming distance between two strings of equal 
-	length is the number of positions at which the corresponding symbols 
+	In information theory, the Hamming distance between two strings of equal
+	length is the number of positions at which the corresponding symbols
 	are different.
 	eg "karolin" and "kathrin" is 3.
 	'''
-	return sum(itertools.imap(operator.ne, str1, str2))
+	return sum(c1 != c2 for c1, c2 in zip(str1, str2))
 
 def intersection(lst1, lst2):
 	lst3 = [value for value in lst1 if value in lst2]
@@ -67,7 +109,7 @@ def get_software_path(tool, abs_path):
 		if(abs_path == ""):
 			sys.exit("ERROR: cannot find "+tool+" in environment; add it to user PATH environment or specify executable using a flag.")
 	if(abs_path != ""):
-		if(os.path.isfile(abs_file)):
+		if(os.path.isfile(abs_path)):
 			tool_path = abs_path
 	return(tool_path)
 
@@ -96,7 +138,10 @@ def check_R_packages(required_packages, R_path):
 	'''
 	Determines whether or not R packages are properly installed
 	'''
-	installed_packages = os.popen(R_path + ''' -e "installed.packages()" | awk '{print $1}' | sort | uniq''').read().strip().split("\n")
+	completed = subprocess.run(
+		[R_path, "--vanilla", "--slave", "-e", "cat(rownames(installed.packages()), sep='\\n')"],
+		stdout=subprocess.PIPE, stderr=subprocess.DEVNULL, text=True)
+	installed_packages = completed.stdout.strip().split("\n")
 	if(not set(required_packages) < set(installed_packages)):
 		sys.exit("ERROR: cannot find the following R package: " + str(set(required_packages) - set(installed_packages)) + "\n" + 
 			"Install it in your R console and then try rerunning proatac (but there may be other missing dependencies).")
@@ -150,7 +195,7 @@ def list_duplicates_of(seq,item):
 
 	# Otherwise figure out .fastq files from directory, do the merging, and infer sample names
 	else:
-		files = os.popen("ls " + input.rstrip("/") +"/*.fastq.gz").read().strip().split("\n")
+		files = sorted(glob.glob(input.rstrip("/") + "/*.fastq.gz"))
 		if(len(files) < 2):
 			sys.exit("ERROR: input determined to be a directory but no paired .fastq.gz files found; QUITTING")
 		files1 = [filename.replace("_R1", "_1").replace("_R2", "_1").replace("_2", "_1") for filename in files]

bap/barcode/cli_barcode.py

@@ -2,22 +2,11 @@
 import os
 import os.path
 import sys
-import shutil
-import yaml
-import random
-import string
-import itertools
-import time
 
-import csv
-import re
-from itertools import groupby
 from ..bapHelp import *
-from pkg_resources import get_distribution
-from subprocess import call, check_call
 
 @click.command()
-@click.version_option()
+@click.version_option(version=bap_version())
 
 @click.argument('mode', type=click.Choice(['v1.0', 'v2.0', 'v2.1', 'v2.1-multi', '10X-v1']))
 
@@ -42,36 +31,37 @@ def main(mode, fastq1, fastq2, fastqi, output, ncores, nreads, nmismatches, reve
 	mode = ['v1.0', 'v2.0', 'v2.1', 'v2.1-multi', '10X-v1'] for bead design\n
 	"""
 
-	__version__ = get_distribution('bap-atac').version
+	__version__ = bap_version()
 	script_dir = os.path.dirname(os.path.realpath(__file__))
 	click.echo(gettime() + "Starting de-barcoding from bap pipeline v%s\n" % __version__)
 
 	# Parse user settings
-	core_call1 = " --fastq1 " + fastq1 + " --fastq2 " + fastq2 + " --ncores " + str(ncores) 
-	core_call2 = " --nreads " + str(nreads) + " --nmismatches " + str(nmismatches) + " --output " + output
-	core_call = core_call1 + core_call2
-
+	core_call = [
+		"--fastq1", fastq1, "--fastq2", fastq2, "--ncores", str(ncores),
+		"--nreads", str(nreads), "--nmismatches", str(nmismatches), "--output", output,
+	]
+
 	# Handle mode to handle the configuration and make the right system call
 	if(mode == "v1.0"):
-		cmd = 'python '+script_dir+'/modes/biorad_v1.py '
-		earlier = " --constant1 " + "TAGCCATCGCATTGC" + " --constant2 " + "TACCTCTGAGCTGAA"
-		later = " --nextera " + "TCGTCGGCAGCGTC" + " --me " + "AGATGTGTATAAGAGACAG"
+		script = script_dir + "/modes/biorad_v1.py"
+		earlier = ["--constant1", "TAGCCATCGCATTGC", "--constant2", "TACCTCTGAGCTGAA"]
+		later = ["--nextera", "TCGTCGGCAGCGTC", "--me", "AGATGTGTATAAGAGACAG"]
 	elif(mode == "v2.0"):
-		cmd = 'python '+script_dir+'/modes/biorad_v2.py '
-		earlier = " --constant1 " + "TATGCATGAC" + " --constant2 " + "AGTCACTGAG"
-		later = " --nextera " + "TGGTAGAGAGGGTG" + " --me " + "AGATGTGTATAAGAGACAG"
+		script = script_dir + "/modes/biorad_v2.py"
+		earlier = ["--constant1", "TATGCATGAC", "--constant2", "AGTCACTGAG"]
+		later = ["--nextera", "TGGTAGAGAGGGTG", "--me", "AGATGTGTATAAGAGACAG"]
 	elif(mode == "v2.1"):
-		cmd = 'python '+script_dir+'/modes/biorad_v2.py '
-		earlier = " --constant1 " + "TATGCATGAC" + " --constant2 " + "AGTCACTGAG"
-		later = " --nextera " + "TCGTCGGCAGCGTC" + " --me " + "AGATGTGTATAAGAGACAG"
+		script = script_dir + "/modes/biorad_v2.py"
+		earlier = ["--constant1", "TATGCATGAC", "--constant2", "AGTCACTGAG"]
+		later = ["--nextera", "TCGTCGGCAGCGTC", "--me", "AGATGTGTATAAGAGACAG"]
 	elif(mode == "v2.1-multi"):
-		cmd = 'python '+script_dir+'/modes/biorad_v2-multi.py '
-		earlier = " --constant1 " + "TATGCATGAC" + " --constant2 " + "AGTCACTGAG"
-		later = " --nextera " + "TCGTCGGCAGCGTC" + " --me " + "AGATGTGTATAAGAGACAG"
+		script = script_dir + "/modes/biorad_v2-multi.py"
+		earlier = ["--constant1", "TATGCATGAC", "--constant2", "AGTCACTGAG"]
+		later = ["--nextera", "TCGTCGGCAGCGTC", "--me", "AGATGTGTATAAGAGACAG"]
 	else:
 		sys.exit(gettime() + "User-supplied mode %s not found!" % mode)
-	
-	# Assemble the final call
-	sys_call = cmd + earlier + later + core_call + " 2> "+output+".stderr.txt"
-	os.system(sys_call)
+
+	# Assemble the final call; stderr is captured to a log file (stdout stays on the terminal)
+	sys_call = [sys.executable, script] + earlier + later + core_call
+	run_cmd(sys_call, log_file=output + ".stderr.txt", log_stderr_only=True)
 

bap/barcode/cli_scaleatac.py

@@ -2,30 +2,19 @@
 import os
 import os.path
 import sys
-import shutil
-import yaml
-import random
-import string
-import itertools
-import time
 import glob
+import gzip
 
-import csv
-import re
-from itertools import groupby
 from ..bapHelp import *
 from .scaleHelp import *
 
-from pkg_resources import get_distribution
-from subprocess import call, check_call
-
 from multiprocessing import Pool, freeze_support
 from Bio import SeqIO
 from Bio.Seq import Seq
 from Bio.SeqIO.QualityIO import FastqGeneralIterator
 
 @click.command()
-@click.version_option()
+@click.version_option(version=bap_version())
 
 @click.option('--fastqs', '-f', help='Path of folder created by mkfastq or bcl2fastq; can be comma separated that will be collapsed into one output.')
 @click.option('--sample', '-s', help='Prefix of the filenames of FASTQs to select; can be comma separated that will be collapsed into one output.')
@@ -43,7 +32,7 @@ def main(fastqs, sample, output, ncores, nreads):
 	Trims, processes Tn5 barcode, and corrects bead barcode in one shot \n
 	"""
 
-	__version__ = get_distribution('bap-atac').version
+	__version__ = bap_version()
 	script_dir = os.path.dirname(os.path.realpath(__file__))
 	click.echo(gettime() + "Starting de-barcoding of scale data from bap pipeline v%s" % __version__)