diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
new file mode 100644
index 0000000..de4feb3
--- /dev/null
+++ b/.github/workflows/CI.yml
@@ -0,0 +1,29 @@
+name: CI
+
+on:
+  push:
+    branches: [master]
+    tags: ['*']
+  pull_request:
+
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version: ['1.10', '1.12']
+        os: [ubuntu-latest, macos-latest]
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.version }}
+      - uses: julia-actions/cache@v2
+      # Manifest.toml is generated with Julia 1.12, so older versions resolve their own within the [compat] bounds
+      - name: Remove Manifest.toml for older Julia versions
+        if: matrix.version != '1.12'
+        run: rm Manifest.toml
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9c62460..3bae9c8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,16 @@
 # CHANGELOG
 
+## v1.10.0
+
+* BUGFIX - `--exclude_regions_file` crashed with "Cannot `convert` an object of type String to an object of type Char" unless `--stranded` was "yes" or "reverse". Added regression tests for both stranded and unstranded runs.
+* FADU is now a proper Julia package (module `FADU`, code in src/, Project.toml and Manifest.toml at the repository root). `julia fadu.jl <options>` still works the same way.
+* FADU's own code is now precompiled (using PrecompileTools and the tiny test dataset as the workload), so runs start about 3x faster than v1.9.2.
+* Adding `Pkg.test()` support (test/runtests.jl) and GitHub Actions CI on Julia v1.10 and v1.12.
+* Adding an experimental `fadu` command for Julia v1.12+ via `Pkg.Apps.add`, and support for `julia -m FADU <options>`.
+* Renamed FADU's internal BAM overlap iterator to `eachalignmentoverlap`. It previously overwrote XAM's `GenomicFeatures.eachoverlap` method, which is not allowed in a precompiled package. Output is unchanged.
+* Moved the tiny regression dataset from test_data/tiny to test/data and the test script to test/run_test.sh.
+* Docker: the container's working directory is now /data (mount your files there) instead of /opt/FADU.
+
 ## v1.9.2
 
 * Raising minimum Julia version to v1.10 and regenerating fadu_pkgs/Manifest.toml under Julia v1.12. Dockerfile now uses the julia:1.12 image.
diff --git a/Dockerfile b/Dockerfile
index 2740f28..6f7391e 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -2,22 +2,19 @@
 FROM julia:1.12
 LABEL maintainer="Shaun Adkins (sadkins@som.umaryland.edu)"
 
-RUN mkdir -p /opt/FADU
-
 WORKDIR /opt/FADU
 
-COPY ./alignment_overlaps.jl /opt/FADU
-COPY ./bam_record.jl /opt/FADU/
-COPY ./fadu.jl /opt/FADU/
-COPY ./feature_counts.jl /opt/FADU/
-COPY ./gff_feature.jl /opt/FADU/
+COPY ./Project.toml ./Manifest.toml ./fadu.jl /opt/FADU/
+COPY ./src /opt/FADU/src
+# The tiny test dataset is used as the precompilation workload
+COPY ./test /opt/FADU/test
+
+# Install and precompile FADU and its dependencies into the image
+ENV JULIA_DEPOT_PATH=/opt/julia_depot
+ENV JULIA_PROJECT=/opt/FADU
+RUN julia -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()'
 
-# force compilation of external, shared packages copied in from Manifest.toml and Project.toml
-# For some reason these would not copy to the "root" home directory
-RUN mkdir -p /tmp/.julia/environments/v1.12
-COPY ./fadu_pkgs/Manifest.toml /tmp/.julia/environments/v1.12/Manifest.toml
-COPY ./fadu_pkgs/Project.toml /tmp/.julia/environments/v1.12/Project.toml
-ENV JULIA_DEPOT_PATH=/tmp/.julia
-RUN julia -e 'using Pkg; Pkg.precompile(); Pkg.instantiate()' && julia -e 'using BGZFStreams, GenomicFeatures, GFF3, Indexes, XAM, BED'
+# Input and output files are expected to be mounted here
+WORKDIR /data
 
-ENTRYPOINT ["julia", "/opt/FADU/fadu.jl"]
\ No newline at end of file
+ENTRYPOINT ["julia", "/opt/FADU/fadu.jl"]
diff --git a/fadu_pkgs/Manifest.toml b/Manifest.toml
similarity index 96%
rename from fadu_pkgs/Manifest.toml
rename to Manifest.toml
index 899d1dd..2d10f1e 100644
--- a/fadu_pkgs/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.12.6"
 manifest_format = "2.0"
-project_hash = "3b4656c2c8a02891df189f30ffd5dbd5b2778bb9"
+project_hash = "314fac2dcffdc3bf9250a8a3c9d875d98c3dacc5"
 
 [[deps.ArgParse]]
 deps = ["Logging", "TextWrap"]
@@ -123,6 +123,12 @@ deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 version = "1.11.0"
 
+[[deps.FADU]]
+deps = ["ArgParse", "BED", "BGZFStreams", "GFF3", "GenomicFeatures", "Indexes", "Logging", "PrecompileTools", "Printf", "StructArrays", "XAM"]
+path = "."
+uuid = "866296ad-388f-4871-8691-d7485d094d7e"
+version = "1.10.0"
+
 [[deps.FASTX]]
 deps = ["Automa", "BioGenerics", "PrecompileTools", "ScanByte", "StringViews", "TranscodingStreams"]
 git-tree-sha1 = "5c4b85ab007cb55d38fc8249ddfab6bf2f48cf06"
@@ -207,9 +213,9 @@ version = "1.8.1"
 
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "07a921781cab75691315adc645096ed5e370cb77"
+git-tree-sha1 = "edbeefc7a4889f528644251bdb5fc9ab5348bc2c"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.3.3"
+version = "1.3.4"
 
 [[deps.Preferences]]
 deps = ["TOML"]
diff --git a/NOTES.md b/NOTES.md
new file mode 100644
index 0000000..d941ae5
--- /dev/null
+++ b/NOTES.md
@@ -0,0 +1,15 @@
+# Notes for package author
+
+## Updating the Manifest.toml file
+
+Manifest.toml pins the dependency versions used for releases, Docker and bioconda. Regenerate it from the repository root with the Julia version you want to target (currently v1.12):
+
+1. `rm Manifest.toml` - This needs to be removed so the next step will resolve a new one
+2. `julia --project=. -e 'using Pkg; Pkg.resolve(); Pkg.instantiate(); Pkg.precompile()'`
+3. Run `julia --project=. -e 'using Pkg; Pkg.test()'` to make sure the output has not changed
+
+## Making a release
+
+1. Bump `version` in Project.toml and `VERSION_NUMBER` in src/FADU.jl (a test fails if they differ)
+2. Add a CHANGELOG.md entry and update "Current FADU release" in README.md and docs/index.md
+3. Commit, then tag with `git tag -a vX.Y.Z -m "vX.Y.Z"` and push the tag
diff --git a/fadu_pkgs/Project.toml b/Project.toml
similarity index 74%
rename from fadu_pkgs/Project.toml
rename to Project.toml
index 2c27ff7..a7e61cf 100644
--- a/fadu_pkgs/Project.toml
+++ b/Project.toml
@@ -1,3 +1,6 @@
+name = "FADU"
+uuid = "866296ad-388f-4871-8691-d7485d094d7e"
+version = "1.10.0"
 authors = ["adkinsrs <sadkins@som.umaryland.edu>"]
 
 [deps]
@@ -8,6 +11,7 @@ GFF3 = "af1dc308-cb6b-11e8-32f0-31192efa90f6"
 GenomicFeatures = "899a7d2d-5c61-547b-bef9-6698a8d05446"
 Indexes = "4ffb77ac-cb80-11e8-1b35-4b78cc642f6d"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 XAM = "d759349c-bcba-11e9-07c2-5b90f8f05f7c"
@@ -19,6 +23,15 @@ BGZFStreams = "0.3.2"
 GFF3 = "0.2.3"
 GenomicFeatures = "2.1.0"
 Indexes = "0.1.3"
+PrecompileTools = "1"
 StructArrays = "0.6.17"
 XAM = "0.3.1"
 julia = "1.10"
+
+[apps.fadu]
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
diff --git a/README.md b/README.md
index 5d1b021..d8143d3 100644
--- a/README.md
+++ b/README.md
@@ -17,15 +17,15 @@ Most current available quantification tools for transcriptomics analyses have be
 
 ## Current FADU release
 
-v1.9.2
+v1.10.0
 
 ## Requirements
 
-Julia - v1.10 or later (fadu_pkgs/Manifest.toml was generated with Julia v1.12)
+Julia - v1.10 or later (Manifest.toml was generated with Julia v1.12)
 
 ### Packages
 
-NOTE: Packages installation instructions can be found at https://docs.julialang.org/en/v1.4/stdlib/Pkg/index.html#Getting-Started-1
+FADU is a Julia package, so these are installed for you by `Pkg.instantiate()` (see "Installation instructions" below).
 
 * ArgParse - v1
 * Logging
@@ -36,6 +36,7 @@ NOTE: Packages installation instructions can be found at https://docs.julialang.
 * StructArrays.jl - v0.6.17
 * XAM.jl - v0.3.1
 * BED.jl - v0.3.0
+* PrecompileTools.jl - v1
 
 ### NOTE on versions
 
@@ -76,8 +77,21 @@ First `cd` to the directory of choice to install FADU
 ```bash
 git clone https://github.com/IGS/FADU.git
 cd FADU
+julia --project=. -e 'using Pkg; Pkg.instantiate()'
 ```
 
+The last step downloads FADU's dependencies and precompiles FADU, which takes a minute or two the first time.  After that, run FADU with `julia ./fadu.jl <options>` from anywhere (use the full path to fadu.jl if you are not in the FADU directory).
+
+### Installing a `fadu` command (Julia v1.12+, experimental)
+
+Julia v1.12 can install FADU as a command-line app:
+
+```bash
+julia -e 'using Pkg; Pkg.Apps.add(url="https://github.com/IGS/FADU.git")'
+```
+
+This puts a `fadu` executable in `~/.julia/bin` (add that directory to your PATH), which takes the same options as `julia ./fadu.jl`.  Julia's app support is still experimental, so the clone-based install above is the recommended route.
+
 ## Example command
 
 For this example we will create a directory called "fadu\_output" to write output into.  This example uses a BAM alignment file of subsampled Wolbachia reads.  This example should take about 60 seconds to run.
@@ -106,13 +120,19 @@ cds1005 816 0.0 0.00    0.00
 
 ## Testing
 
-A small regression dataset (a ~20 kb slice of an E. coli paired-end RNA-Seq BAM, plus the matching GFF3 annotation) lives in "test_data/tiny".  To run FADU on it and compare against the expected output:
+A small regression dataset (a ~20 kb slice of an E. coli paired-end RNA-Seq BAM, plus the matching GFF3 annotation) lives in "test/data".  To run the test suite:
+
+```bash
+julia --project=. -e 'using Pkg; Pkg.test()'
+```
+
+To check only the command-line output against the expected counts, without Julia's test framework:
 
 ```bash
-bash ./test_data/tiny/run_test.sh
+bash ./test/run_test.sh
 ```
 
-By default the script runs `julia --project=./fadu_pkgs ./fadu.jl`.  Set the `FADU_CMD` environment variable to test a different FADU command (e.g. an installed copy).
+By default the script runs `julia ./fadu.jl`.  Set the `FADU_CMD` environment variable to test a different FADU command (e.g. an installed copy).
 
 ## Usage
 
@@ -218,10 +238,10 @@ docker build -t fadu .
 This Docker image uses Julia 1.12 and installs the various packages required to run FADU. Once the image has been built you can create and run a FADU container by running the following:
 
 ```bash
-docker run -v $PWD:/opt/FADU fadu <fadu_options>
+docker run -v $PWD:/data fadu <fadu_options>
 ```
 
-This command assumes that your input files are in your current working directory and that your output will be written to a space relative to your current working directory as well.
+This command mounts your current working directory at /data, which is the container's working directory, so input and output paths relative to your current working directory work as-is.
 
 ## Issues
 
diff --git a/docs/index.md b/docs/index.md
index 85cbeee..5e6450a 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -17,15 +17,15 @@ Most current available quantification tools for transcriptomics analyses have be
 
 ## Current FADU release
 
-v1.9.2
+v1.10.0
 
 ## Requirements
 
-Julia - v1.10 or later (fadu_pkgs/Manifest.toml was generated with Julia v1.12)
+Julia - v1.10 or later (Manifest.toml was generated with Julia v1.12)
 
 ### Packages
 
-NOTE: Packages installation instructions can be found at https://docs.julialang.org/en/v1.4/stdlib/Pkg/index.html#Getting-Started-1
+FADU is a Julia package, so these are installed for you by `Pkg.instantiate()` (see "Installation instructions" below).
 
 * ArgParse - v1
 * Logging
@@ -36,6 +36,7 @@ NOTE: Packages installation instructions can be found at https://docs.julialang.
 * StructArrays.jl - v0.6.17
 * XAM.jl - v0.3.1
 * BED.jl - v0.3.0
+* PrecompileTools.jl - v1
 
 ### NOTE on versions
 
@@ -76,8 +77,21 @@ First `cd` to the directory of choice to install FADU
 ```bash
 git clone https://github.com/IGS/FADU.git
 cd FADU
+julia --project=. -e 'using Pkg; Pkg.instantiate()'
 ```
 
+The last step downloads FADU's dependencies and precompiles FADU, which takes a minute or two the first time.  After that, run FADU with `julia ./fadu.jl <options>` from anywhere (use the full path to fadu.jl if you are not in the FADU directory).
+
+### Installing a `fadu` command (Julia v1.12+, experimental)
+
+Julia v1.12 can install FADU as a command-line app:
+
+```bash
+julia -e 'using Pkg; Pkg.Apps.add(url="https://github.com/IGS/FADU.git")'
+```
+
+This puts a `fadu` executable in `~/.julia/bin` (add that directory to your PATH), which takes the same options as `julia ./fadu.jl`.  Julia's app support is still experimental, so the clone-based install above is the recommended route.
+
 ## Example command
 
 For this example we will create a directory called "fadu\_output" to write output into.  This example uses a BAM alignment file of subsampled Wolbachia reads. This example should take about 60 seconds to run.
@@ -107,13 +121,19 @@ cds1005 816 0.0 0.00    0.00
 
 ## Testing
 
-A small regression dataset (a ~20 kb slice of an E. coli paired-end RNA-Seq BAM, plus the matching GFF3 annotation) lives in "test_data/tiny".  To run FADU on it and compare against the expected output:
+A small regression dataset (a ~20 kb slice of an E. coli paired-end RNA-Seq BAM, plus the matching GFF3 annotation) lives in "test/data".  To run the test suite:
+
+```bash
+julia --project=. -e 'using Pkg; Pkg.test()'
+```
+
+To check only the command-line output against the expected counts, without Julia's test framework:
 
 ```bash
-bash ./test_data/tiny/run_test.sh
+bash ./test/run_test.sh
 ```
 
-By default the script runs `julia --project=./fadu_pkgs ./fadu.jl`.  Set the `FADU_CMD` environment variable to test a different FADU command (e.g. an installed copy).
+By default the script runs `julia ./fadu.jl`.  Set the `FADU_CMD` environment variable to test a different FADU command (e.g. an installed copy).
 
 ## Usage
 
@@ -219,10 +239,10 @@ docker build -t fadu .
 This Docker image uses Julia 1.12 and installs the various packages required to run FADU. Once the image has been built you can create and run a FADU container by running the following:
 
 ```bash
-docker run -v $PWD:/opt/FADU fadu <fadu_options>
+docker run -v $PWD:/data fadu <fadu_options>
 ```
 
-This command assumes that your input files are in your current working directory and that your output will be written to a space relative to your current working directory as well.
+This command mounts your current working directory at /data, which is the container's working directory, so input and output paths relative to your current working directory work as-is.
 
 ## Issues
 
diff --git a/fadu.jl b/fadu.jl
old mode 100644
new mode 100755
index 50c4380..d84ad9e
--- a/fadu.jl
+++ b/fadu.jl
@@ -1,304 +1,7 @@
 #!/usr/bin/env julia
 
-"""
-FADU.jl - Feature Aggregate Depth Utility.
-
-Description - Generate fractional counts of alignments that map to non-overlapping portions of genes
-Requires (version number listed is earliest supported version):
-    Julia - v1.10
-    * Packages *
-    BGZFStreams.jl - v0.3.2
-    GenomicFeatures.jl - v2.1.0
-    GFF3 - v0.2.3
-    Indexes.jl - v0.1.3
-    StructArrays.jl - v0.6.17
-    XAM.jl - v0.3.1
-    BED.jl - v0.3.0
-
-By: Shaun Adkins (sadkins@som.umaryland.edu)
-    Matthew Chung (mattchung@umaryland.edu)
-
-"""
-
-# The macro on modules and functions makes the code available to all worker processes
-using ArgParse
-using Logging
-using XAM: BAM, SAM
-using GenomicFeatures
-using GFF3
-using Indexes
-using Printf
-using StructArrays
-using BED
-
-include("alignment_overlaps.jl")
-include("bam_record.jl")
-include("feature_counts.jl")
-include("gff_feature.jl")
-
-const VERSION_NUMBER = "1.9.2"    # Version number of the FADU program
-const MAX_FRAGMENT_SIZE = 1000 # Maximum size of fragment.  If exceeded, fragment will be considered two reads
-const EM_ITER_DEFAULT = 1 # Number of iterations to do EM-algorithm
-
-### Functions to clean up, then move to their proper "include" script
-
-function remove_excluded_regions!(uniq_coords::Dict{String, Dict}, excluded_regions_file, stranded::Bool=false)
-    """Remove any coordinates that overlap the regions in the BED file. Edits "uniq_coords" in-place."""
-
-    @info("Removing alignments that overlap excluded regions...")
-    # if file does not exist, throw error
-    isfile(excluded_regions_file) || throw(SystemError("Excluded regions file does not seem to exist. Please check supplied path."))
-
-    exclude_regions = open(collect, BED.Reader, excluded_regions_file)
-    exclude_coords = Dict{String, Dict}()
-    # Read in all excluded regions, by strand if the "stranded" argument is passed
-    for region in exclude_regions
-        # Store chrom since that should match seqid
-        seqid = BED.chrom(region)
-        if !haskey(uniq_coords, seqid)
-            continue
-        end
-
-        # "get!" will create the key if it doesn't exist in addition to returning the value
-        seqid_exclude_coords = get!(exclude_coords, seqid, Dict{Char, BitSet}())
-
-        start = BED.chromstart(region)
-        stop = BED.chromend(region)
-
-        # Store coordinates to exclude
-        # If unstranded, then everything goes in "+"
-        if stranded
-            strand = convert(Char, BED.strand(region))
-            union!(get!(seqid_exclude_coords, strand, BitSet()), BitSet(start:stop))
-        else
-            union!(get!(seqid_exclude_coords, "+", BitSet()), BitSet(start:stop))
-        end
-    end
-
-
-    # Remove any alignments that overlap these regions
-    for (seqid, strand_coords) in exclude_coords
-        if haskey(uniq_coords, seqid)
-            for (strand, coords) in strand_coords
-                if haskey(uniq_coords[seqid], strand)
-                    setdiff!(uniq_coords[seqid][strand], coords)
-                end
-            end
-        end
-    end
-
-end
-
-
-########
-# Main #
-########
-
-function parse_commandline()
-    s = ArgParseSettings(description = "Generate counts of reads that map to non-overlapping portions of genes",
-    prog = "fadu.jl",
-	version = VERSION_NUMBER,
-	add_version = true,
-    add_help = true)
-
-    # The macro to add a table of arguments and options to the ArgParseSettings object
-    @add_arg_table! s begin
-        "--bam_file", "-b"
-            help = "Path to BAM file (SAM is not supported)."
-            metavar = "/path/to/file.bam"
-            required = true
-        "--gff3_file", "-g"
-            help = "Path to GFF3-formatted annotation file (GTF is not supported)."
-            metavar = "/path/to/annotation.gff3"
-            required = true
-        "--output_dir", "-o"
-            help = "Directory to write the output."
-            metavar = "/path/to/output/dir/"
-            required = true
-        "--no_output_header"
-            help = "If enabled, do not write the header line to the output file."
-            action = :store_true
-            dest_name = "no_header"
-        "--exclude_regions_file", "-x"
-            help = "Path to BED file containing regions to exclude from analysis.  Any alignments that overlap these regions will be ignored."
-            metavar = "/path/to/regions.bed"
-        "--stranded", "-s"
-            help = "Indicate if BAM reads are from a strand-specific assay. Choose between 'yes', 'no', or 'reverse'."
-            default = "no"
-            range_tester = (x->x in ["yes", "no", "reverse"])
-        "--feature_type", "-f"
-            help = "Which GFF3/GTF feature type (column 3) to obtain. readcount statistics for. Case-sensitive."
-            default = "gene"
-        "--attribute_type", "-a"
-            help = "Which GFF3/GTF feature type (column 9) to obtain. readcount statistics for. Case-sensitive."
-            default = "ID"
-        "--keep_only_proper_pairs", "-p"
-            help = "If enabled, keep only properly paired reads when performing calculations."
-            action = :store_true
-            dest_name = "pp_only"
-        "--max_fragment_size", "-m"
-            help = "If the fragment size of properly-paired reads exceeds this value, process pair as single reads instead of as a fragment. Setting this value to 0 will make every fragment pair be processed as two individual reads. If --keep_only_proper_pairs is enabled, then any fragment exceeding this value will be discarded."
-            default = MAX_FRAGMENT_SIZE
-            arg_type = Int
-            range_tester = (x->typemin(UInt)<=x<=typemax(UInt))
-        "--remove_multimapped", "-M"
-            help = "If enabled, remove any reads or fragments that are mapped to multiple regions of the genome, indiated by the 'NH' attribute being greater than 1.  Otherwise, use EM algorithm to quantify reads after all other reads are counted."
-            action = :store_true
-            dest_name = "rm_multimap"
-        "--em_iterations", "-e"
-            help = "Number of iterations to perform EM algorithm on multimapped read depth. Only applies if --remove_multimapped flag is not passed (is disabled)"
-            default = EM_ITER_DEFAULT
-            arg_type = Int
-            range_tester = (x->x>0)
-            dest_name = "em_iter"
-        "--log_level"
-            help = "Set the log level.  Options are: DEBUG, INFO, WARNING, ERROR, CRITICAL."
-            default = "INFO"
-            metavar = "LEVEL"
-            range_tester = (x->x in ["DEBUG", "INFO", "WARNING", "ERROR"])
-        "--log_file"
-            help = "Path to log file.  If not specified, log messages will be printed to STDERR."
-            metavar = "/path/to/logfile.log"
-
-    # Will not add log_file or debug options for now
-    end
-    # Converts the ArgParseSettings object into key/value pairs
-    return parse_args(s)
-end
-
-function setup_logger(args)
-    """Set up logging."""
-    stream = args["log_file"] === nothing ? stderr : open(args["log_file"], "w")
-    logleveldict = Dict("DEBUG" => -1000, "INFO" => 0, "WARNING" => 1000, "ERROR" => 2000)
-    loglevel = get(logleveldict, uppercase(args["log_level"]), 0)
-    logger = ConsoleLogger(stream, loglevel)
-    global_logger(logger)
-    # TODO: prevent debug messages from printing the file and line number
-
-end
-
-function validate_args(args)
-    """Validate the passed arguments."""
-    @info("Validating arguments...")
-    isfile(args["bam_file"]) || error("BAM file does not exist. Please check supplied path")
-    isfile(args["gff3_file"]) || error("GFF3 file does not exist. Please check supplied path")
-    if !isdir(args["output_dir"])
-        @info("Creating output directory at ", args["output_dir"])
-        mkdir(args["output_dir"])
-    end
-end
-
-function main()
-    args = parse_commandline()
-    setup_logger(args)
-
-    try
-        validate_args(args)
-    catch e
-        @error(e)
-        exit(1)
-    end
-
-    @info("Parsed args:")
-    for (arg,val) in args
-        @info("  $arg  =>  $val")
-    end
-
-    @info("Processing annotation features...")
-    features = try
-        open(collect, GFF3.Reader, args["gff3_file"])    # Array of GFF3.Record objects
-    catch e
-        println(stderr, "ERROR: GFF3 file was unable to be read, due to possible improper formatting.  Please consult https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md for how to properly format the GFF3 file.  Actual error is below.  Program will now exit.\n")
-        showerror(stderr, e)
-        exit(1)
-    end
-
-    # Only keep features for the chosen feature type
-    filter!(x -> GFF3.featuretype(x) == args["feature_type"], features)
-
-    @info("Getting unique coordinates per contig and feature...")
-    uniq_coords = create_uniq_coords_dict(features, args["stranded"])
-
-    # If exclude_regions is passed, then remove any alignments that overlap these regions
-    if args["exclude_regions_file"] !== nothing
-        try
-            remove_excluded_regions!(uniq_coords, args["exclude_regions_file"], isstranded(args["stranded"]))
-        catch e
-            println(stderr, "ERROR: Excluded regions file was unable to be read, due to possible improper formatting.  Please consult https://genome.ucsc.edu/FAQ/FAQformat.html#format1 for how to properly format the BED file.  Actual error is below.  Program will now exit.\n")
-            showerror(stderr, e)
-            exit(1)
-        end
-    end
-
-    @info("Initializing dictionary of feature count information")
-    feat_overlaps = create_feat_overlaps_dict(features, uniq_coords, args["attribute_type"], args["stranded"])
-
-    bai_file = replace(args["bam_file"], ".bam" => ".bai")
-    # Attempt to find index file for BAM file in same directory
-    if !isfile(bai_file)
-        bai_file = args["bam_file"] * ".bai"
-        if !isfile(bai_file)
-            bam = args["bam_file"]
-            error("Attempted to find index file for BAM file $bam but could not find one")
-        end
-    end
-    @debug("Found BAM index file at $bai_file")
-
-    # Open a BAM file and iterate over records overlapping GFF features.
-    @info("Opening BAM alignment file...")
-    reader = open(BAM.Reader, args["bam_file"], index = bai_file)
-
-    @info("Now finding overlaps between alignment and annotation records...")
-    multimapped_dict = process_all_feature_overlaps(feat_overlaps, features, reader, args)
-    close(reader)
-
-    # If multimapped reads are kept, use EM algorithm to count and re-add back into the feature counts
-    if !args["rm_multimap"]
-        num_multimaps = length(collect(keys(multimapped_dict)))
-        @debug("Multimapped alignment templates: ", num_multimaps)
-
-        featurenames = collect(keys(feat_overlaps))
-        mm_feat_overlaps = Dict{String, FeatureOverlap}(featurename => initialize_overlap_info(coordinate_set(feat_overlaps[featurename])) for featurename in featurenames)
-
-        @info("Counting and adjusting multimapped alignment feature counts via Expectation-Maximization algorithm...")
-        while args["em_iter"] > 0
-            @debug("\tEM iterations left: ", args["em_iter"])
-            adjusted_overlaps = merge_mm_counts(feat_overlaps, mm_feat_overlaps, false)
-            args["em_iter"] -= 1
-            mm_feat_overlaps = compute_mm_counts_by_em(adjusted_overlaps, multimapped_dict)
-        end
-        # Last iteration the alignment counts are added too
-        merge_mm_counts!(feat_overlaps, mm_feat_overlaps, true)
-    end
-
-    # Calculate sum of all counts for each feature
-    totalcounts = calc_totalcounts(feat_overlaps)
-
-    @info("Writing counts output to file...")
-    out_file = joinpath(args["output_dir"], splitext(basename(args["bam_file"]))[1]) * ".counts.txt"
-    out_f = open(out_file, "w")
-    if !args["no_header"]
-        write(out_f, "featureID\tuniq_len\tnum_alignments\tcounts\ttpm\n")
-    end
-    # Write output, sorted alphabetically
-    for featurename in sort(collect(keys(feat_overlaps)))
-        uniq_len::UInt = length(coordinate_set(feat_overlaps[featurename]))
-        num_alignments = totalalignments(feat_overlaps[featurename])
-        feat_counts = featurecounts(feat_overlaps[featurename])
-        if isnan(feat_counts)
-            feat_counts = zero(Float32)
-        end
-        tpm = calc_tpm(uniq_len, totalcounts, feat_counts)
-        if isnan(tpm)
-            tpm = zero(Float32)
-        end
-        s = @sprintf("%s\t%i\t%.1f\t%.2f\t%.2f\n", featurename, uniq_len, num_alignments, feat_counts, tpm)
-        write(out_f, s)
-    end
-    close(out_f)
-
-    @info("FADU is complete!  Exiting!")
-end
-
-main()
-exit()
+# Launcher for FADU. Loads the FADU package from this directory, so the
+# command `julia fadu.jl [options]` keeps working from a clone of the repository.
+pushfirst!(LOAD_PATH, @__DIR__)
+using FADU
+exit(FADU.main(ARGS))
diff --git a/fadu_pkgs/NOTES.md b/fadu_pkgs/NOTES.md
deleted file mode 100644
index 8f0a111..0000000
--- a/fadu_pkgs/NOTES.md
+++ /dev/null
@@ -1,9 +0,0 @@
-# Notes for package author
-
-## Updating the Manifest.toml file
-
-Project.toml is a plain Julia environment (no name/uuid/version), so the Manifest can be regenerated directly from the repository root with the Julia version you want to target (currently v1.12):
-
-1. `rm fadu_pkgs/Manifest.toml` - This needs to be removed so the next step will resolve a new one
-2. `julia --project=fadu_pkgs -e 'using Pkg; Pkg.resolve(); Pkg.instantiate(); Pkg.precompile()'`
-3. Run `bash test_data/tiny/run_test.sh` to make sure the output has not changed
diff --git a/src/FADU.jl b/src/FADU.jl
new file mode 100644
index 0000000..e20f0d7
--- /dev/null
+++ b/src/FADU.jl
@@ -0,0 +1,334 @@
+"""
+FADU.jl - Feature Aggregate Depth Utility.
+
+Description - Generate fractional counts of alignments that map to non-overlapping portions of genes
+Requires (version number listed is earliest supported version):
+    Julia - v1.10
+    * Packages *
+    ArgParse.jl - v1
+    BGZFStreams.jl - v0.3.2
+    GenomicFeatures.jl - v2.1.0
+    GFF3 - v0.2.3
+    Indexes.jl - v0.1.3
+    StructArrays.jl - v0.6.17
+    XAM.jl - v0.3.1
+    BED.jl - v0.3.0
+    PrecompileTools.jl - v1
+
+By: Shaun Adkins (sadkins@som.umaryland.edu)
+    Matthew Chung (mattchung@umaryland.edu)
+
+"""
+module FADU
+
+# `main` is deliberately not exported: it is marked with `@main`, so exporting it would
+# make any script that runs `using FADU` also run FADU on exit. Call it as `FADU.main`.
+
+using ArgParse
+using Logging
+using XAM: BAM, SAM
+using GenomicFeatures
+using GFF3
+using Indexes
+using Printf
+using StructArrays
+using BED
+using PrecompileTools
+
+include("alignment_overlaps.jl")
+include("bam_record.jl")
+include("feature_counts.jl")
+include("gff_feature.jl")
+
+const VERSION_NUMBER = "1.10.0"    # Version number of the FADU program. Keep in sync with Project.toml
+const MAX_FRAGMENT_SIZE = 1000 # Maximum size of fragment.  If exceeded, fragment will be considered two reads
+const EM_ITER_DEFAULT = 1 # Number of iterations to do EM-algorithm
+
+### Functions to clean up, then move to their proper "include" script
+
+function remove_excluded_regions!(uniq_coords::Dict{String, Dict}, excluded_regions_file, stranded::Bool=false)
+    """Remove any coordinates that overlap the regions in the BED file. Edits "uniq_coords" in-place."""
+
+    @info("Removing alignments that overlap excluded regions...")
+    # if file does not exist, throw error
+    isfile(excluded_regions_file) || throw(SystemError("Excluded regions file does not seem to exist. Please check supplied path."))
+
+    exclude_regions = open(collect, BED.Reader, excluded_regions_file)
+    exclude_coords = Dict{String, Dict}()
+    # Read in all excluded regions, by strand if the "stranded" argument is passed
+    for region in exclude_regions
+        # Store chrom since that should match seqid
+        seqid = BED.chrom(region)
+        if !haskey(uniq_coords, seqid)
+            continue
+        end
+
+        # "get!" will create the key if it doesn't exist in addition to returning the value
+        seqid_exclude_coords = get!(exclude_coords, seqid, Dict{Char, BitSet}())
+
+        start = BED.chromstart(region)
+        stop = BED.chromend(region)
+
+        # Store coordinates to exclude
+        # If unstranded, then everything goes in "+"
+        if stranded
+            strand = convert(Char, BED.strand(region))
+            union!(get!(seqid_exclude_coords, strand, BitSet()), BitSet(start:stop))
+        else
+            union!(get!(seqid_exclude_coords, '+', BitSet()), BitSet(start:stop))
+        end
+    end
+
+
+    # Remove any alignments that overlap these regions
+    for (seqid, strand_coords) in exclude_coords
+        if haskey(uniq_coords, seqid)
+            for (strand, coords) in strand_coords
+                if haskey(uniq_coords[seqid], strand)
+                    setdiff!(uniq_coords[seqid][strand], coords)
+                end
+            end
+        end
+    end
+
+end
+
+
+########
+# Main #
+########
+
+function parse_commandline(cli_args::AbstractVector{<:AbstractString})
+    s = ArgParseSettings(description = "Generate counts of reads that map to non-overlapping portions of genes",
+    prog = "fadu.jl",
+	version = VERSION_NUMBER,
+	add_version = true,
+    add_help = true)
+
+    # The macro to add a table of arguments and options to the ArgParseSettings object
+    @add_arg_table! s begin
+        "--bam_file", "-b"
+            help = "Path to BAM file (SAM is not supported)."
+            metavar = "/path/to/file.bam"
+            required = true
+        "--gff3_file", "-g"
+            help = "Path to GFF3-formatted annotation file (GTF is not supported)."
+            metavar = "/path/to/annotation.gff3"
+            required = true
+        "--output_dir", "-o"
+            help = "Directory to write the output."
+            metavar = "/path/to/output/dir/"
+            required = true
+        "--no_output_header"
+            help = "If enabled, do not write the header line to the output file."
+            action = :store_true
+            dest_name = "no_header"
+        "--exclude_regions_file", "-x"
+            help = "Path to BED file containing regions to exclude from analysis.  Any alignments that overlap these regions will be ignored."
+            metavar = "/path/to/regions.bed"
+        "--stranded", "-s"
+            help = "Indicate if BAM reads are from a strand-specific assay. Choose between 'yes', 'no', or 'reverse'."
+            default = "no"
+            range_tester = (x->x in ["yes", "no", "reverse"])
+        "--feature_type", "-f"
+            help = "Which GFF3/GTF feature type (column 3) to obtain. readcount statistics for. Case-sensitive."
+            default = "gene"
+        "--attribute_type", "-a"
+            help = "Which GFF3/GTF feature type (column 9) to obtain. readcount statistics for. Case-sensitive."
+            default = "ID"
+        "--keep_only_proper_pairs", "-p"
+            help = "If enabled, keep only properly paired reads when performing calculations."
+            action = :store_true
+            dest_name = "pp_only"
+        "--max_fragment_size", "-m"
+            help = "If the fragment size of properly-paired reads exceeds this value, process pair as single reads instead of as a fragment. Setting this value to 0 will make every fragment pair be processed as two individual reads. If --keep_only_proper_pairs is enabled, then any fragment exceeding this value will be discarded."
+            default = MAX_FRAGMENT_SIZE
+            arg_type = Int
+            range_tester = (x->typemin(UInt)<=x<=typemax(UInt))
+        "--remove_multimapped", "-M"
+            help = "If enabled, remove any reads or fragments that are mapped to multiple regions of the genome, indiated by the 'NH' attribute being greater than 1.  Otherwise, use EM algorithm to quantify reads after all other reads are counted."
+            action = :store_true
+            dest_name = "rm_multimap"
+        "--em_iterations", "-e"
+            help = "Number of iterations to perform EM algorithm on multimapped read depth. Only applies if --remove_multimapped flag is not passed (is disabled)"
+            default = EM_ITER_DEFAULT
+            arg_type = Int
+            range_tester = (x->x>0)
+            dest_name = "em_iter"
+        "--log_level"
+            help = "Set the log level.  Options are: DEBUG, INFO, WARNING, ERROR, CRITICAL."
+            default = "INFO"
+            metavar = "LEVEL"
+            range_tester = (x->x in ["DEBUG", "INFO", "WARNING", "ERROR"])
+        "--log_file"
+            help = "Path to log file.  If not specified, log messages will be printed to STDERR."
+            metavar = "/path/to/logfile.log"
+
+    # Will not add log_file or debug options for now
+    end
+    # Converts the ArgParseSettings object into key/value pairs
+    return parse_args(cli_args, s)
+end
+
+function setup_logger(args)
+    """Set up logging."""
+    stream = args["log_file"] === nothing ? stderr : open(args["log_file"], "w")
+    logleveldict = Dict("DEBUG" => -1000, "INFO" => 0, "WARNING" => 1000, "ERROR" => 2000)
+    loglevel = get(logleveldict, uppercase(args["log_level"]), 0)
+    logger = ConsoleLogger(stream, loglevel)
+    global_logger(logger)
+    # TODO: prevent debug messages from printing the file and line number
+
+end
+
+function validate_args(args)
+    """Validate the passed arguments."""
+    @info("Validating arguments...")
+    isfile(args["bam_file"]) || error("BAM file does not exist. Please check supplied path")
+    isfile(args["gff3_file"]) || error("GFF3 file does not exist. Please check supplied path")
+    if !isdir(args["output_dir"])
+        @info("Creating output directory at ", args["output_dir"])
+        mkdir(args["output_dir"])
+    end
+end
+
+"""
+    main(cli_args=ARGS) -> Int
+
+Run FADU with command-line style arguments and return the exit code.
+"""
+function main(cli_args::AbstractVector{<:AbstractString}=ARGS)
+    args = parse_commandline(cli_args)
+    setup_logger(args)
+
+    try
+        validate_args(args)
+    catch e
+        @error(e)
+        return 1
+    end
+
+    @info("Parsed args:")
+    for (arg,val) in args
+        @info("  $arg  =>  $val")
+    end
+
+    @info("Processing annotation features...")
+    features = try
+        open(collect, GFF3.Reader, args["gff3_file"])    # Array of GFF3.Record objects
+    catch e
+        println(stderr, "ERROR: GFF3 file was unable to be read, due to possible improper formatting.  Please consult https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md for how to properly format the GFF3 file.  Actual error is below.  Program will now exit.\n")
+        showerror(stderr, e)
+        return 1
+    end
+
+    # Only keep features for the chosen feature type
+    filter!(x -> GFF3.featuretype(x) == args["feature_type"], features)
+
+    @info("Getting unique coordinates per contig and feature...")
+    uniq_coords = create_uniq_coords_dict(features, args["stranded"])
+
+    # If exclude_regions is passed, then remove any alignments that overlap these regions
+    if args["exclude_regions_file"] !== nothing
+        try
+            remove_excluded_regions!(uniq_coords, args["exclude_regions_file"], isstranded(args["stranded"]))
+        catch e
+            println(stderr, "ERROR: Excluded regions file was unable to be read, due to possible improper formatting.  Please consult https://genome.ucsc.edu/FAQ/FAQformat.html#format1 for how to properly format the BED file.  Actual error is below.  Program will now exit.\n")
+            showerror(stderr, e)
+            return 1
+        end
+    end
+
+    @info("Initializing dictionary of feature count information")
+    feat_overlaps = create_feat_overlaps_dict(features, uniq_coords, args["attribute_type"], args["stranded"])
+
+    bai_file = replace(args["bam_file"], ".bam" => ".bai")
+    # Attempt to find index file for BAM file in same directory
+    if !isfile(bai_file)
+        bai_file = args["bam_file"] * ".bai"
+        if !isfile(bai_file)
+            bam = args["bam_file"]
+            error("Attempted to find index file for BAM file $bam but could not find one")
+        end
+    end
+    @debug("Found BAM index file at $bai_file")
+
+    # Open a BAM file and iterate over records overlapping GFF features.
+    @info("Opening BAM alignment file...")
+    reader = open(BAM.Reader, args["bam_file"], index = bai_file)
+
+    @info("Now finding overlaps between alignment and annotation records...")
+    multimapped_dict = process_all_feature_overlaps(feat_overlaps, features, reader, args)
+    close(reader)
+
+    # If multimapped reads are kept, use EM algorithm to count and re-add back into the feature counts
+    if !args["rm_multimap"]
+        num_multimaps = length(collect(keys(multimapped_dict)))
+        @debug("Multimapped alignment templates: ", num_multimaps)
+
+        featurenames = collect(keys(feat_overlaps))
+        mm_feat_overlaps = Dict{String, FeatureOverlap}(featurename => initialize_overlap_info(coordinate_set(feat_overlaps[featurename])) for featurename in featurenames)
+
+        @info("Counting and adjusting multimapped alignment feature counts via Expectation-Maximization algorithm...")
+        while args["em_iter"] > 0
+            @debug("\tEM iterations left: ", args["em_iter"])
+            adjusted_overlaps = merge_mm_counts(feat_overlaps, mm_feat_overlaps, false)
+            args["em_iter"] -= 1
+            mm_feat_overlaps = compute_mm_counts_by_em(adjusted_overlaps, multimapped_dict)
+        end
+        # Last iteration the alignment counts are added too
+        merge_mm_counts!(feat_overlaps, mm_feat_overlaps, true)
+    end
+
+    # Calculate sum of all counts for each feature
+    totalcounts = calc_totalcounts(feat_overlaps)
+
+    @info("Writing counts output to file...")
+    out_file = joinpath(args["output_dir"], splitext(basename(args["bam_file"]))[1]) * ".counts.txt"
+    out_f = open(out_file, "w")
+    if !args["no_header"]
+        write(out_f, "featureID\tuniq_len\tnum_alignments\tcounts\ttpm\n")
+    end
+    # Write output, sorted alphabetically
+    for featurename in sort(collect(keys(feat_overlaps)))
+        uniq_len::UInt = length(coordinate_set(feat_overlaps[featurename]))
+        num_alignments = totalalignments(feat_overlaps[featurename])
+        feat_counts = featurecounts(feat_overlaps[featurename])
+        if isnan(feat_counts)
+            feat_counts = zero(Float32)
+        end
+        tpm = calc_tpm(uniq_len, totalcounts, feat_counts)
+        if isnan(tpm)
+            tpm = zero(Float32)
+        end
+        s = @sprintf("%s\t%i\t%.1f\t%.2f\t%.2f\n", featurename, uniq_len, num_alignments, feat_counts, tpm)
+        write(out_f, s)
+    end
+    close(out_f)
+
+    @info("FADU is complete!  Exiting!")
+    return 0
+end
+
+# Mark main as the entry point for `julia -m FADU` and Pkg apps (Julia 1.11+)
+@static if isdefined(Base, Symbol("@main"))
+    @main
+end
+
+# Run FADU on the tiny test dataset during precompilation so the first real run starts fast.
+# Skipped if the test data was not shipped alongside the source.
+@setup_workload begin
+    data_dir = joinpath(pkgdir(FADU), "test", "data")
+    if isfile(joinpath(data_dir, "tiny.bam"))
+        @compile_workload begin
+            with_logger(NullLogger()) do
+                mktempdir() do out_dir
+                    main(["-g", joinpath(data_dir, "tiny.gff3"), "-b", joinpath(data_dir, "tiny.bam"),
+                          "-o", out_dir, "-s", "reverse", "-f", "CDS", "-a", "ID"])
+                end
+            end
+        end
+    end
+end
+
+end # module FADU
diff --git a/alignment_overlaps.jl b/src/alignment_overlaps.jl
similarity index 81%
rename from alignment_overlaps.jl
rename to src/alignment_overlaps.jl
index d84acee..b826d58 100644
--- a/alignment_overlaps.jl
+++ b/src/alignment_overlaps.jl
@@ -12,6 +12,7 @@ By: Shaun Adkins (sadkins@som.umaryland.edu)
 import BGZFStreams
 
 ### OverlapIterator - Extended version of original BAM OverlapIterator
+# Kept separate from GenomicFeatures.eachoverlap so it does not overwrite the XAM method of the same signature
 
 struct OverlapIterator{T}
     reader::BAM.Reader{T}
@@ -21,15 +22,15 @@ struct OverlapIterator{T}
     max_frag_size::UInt
 end
 
-function GenomicFeatures.eachoverlap(reader::BAM.Reader, interval::Interval, strand_type::String="no", max_frag_size::UInt=1000)
-    return GenomicFeatures.eachoverlap(reader, interval.seqname, interval.first:interval.last, strand_type, max_frag_size)
+function eachalignmentoverlap(reader::BAM.Reader, interval::Interval, strand_type::String="no", max_frag_size::UInt=1000)
+    return eachalignmentoverlap(reader, interval.seqname, interval.first:interval.last, strand_type, max_frag_size)
 end
 
-function GenomicFeatures.eachoverlap(reader::BAM.Reader, interval, strand_type::String="no", max_frag_size::UInt=1000)
-    return GenomicFeatures.eachoverlap(reader, convert(Interval, interval), strand_type, max_frag_size)
+function eachalignmentoverlap(reader::BAM.Reader, interval, strand_type::String="no", max_frag_size::UInt=1000)
+    return eachalignmentoverlap(reader, convert(Interval, interval), strand_type, max_frag_size)
 end
 
-function GenomicFeatures.eachoverlap(reader::BAM.Reader, refname::AbstractString, interval::UnitRange, strand_type::String="no", max_frag_size::UInt=1000)
+function eachalignmentoverlap(reader::BAM.Reader, refname::AbstractString, interval::UnitRange, strand_type::String="no", max_frag_size::UInt=1000)
     return OverlapIterator(reader, String(refname), interval, strand_type, max_frag_size)
 end
 
diff --git a/bam_record.jl b/src/bam_record.jl
similarity index 100%
rename from bam_record.jl
rename to src/bam_record.jl
diff --git a/feature_counts.jl b/src/feature_counts.jl
similarity index 99%
rename from feature_counts.jl
rename to src/feature_counts.jl
index eecbdfa..037b980 100644
--- a/feature_counts.jl
+++ b/src/feature_counts.jl
@@ -177,7 +177,7 @@ function process_feature_overlaps!(feat_overlaps::Dict{String, FeatureOverlap},
     count = 0
 
     # Returning SuperBAMRecord objects instead of BAM.Record objects
-    for rec in eachoverlap(reader, feature, args["stranded"], max_frag_size)
+    for rec in eachalignmentoverlap(reader, feature, args["stranded"], max_frag_size)
         alignmentinterval = interval(rec)
         @assert alignmentinterval !== nothing "All iterated records should be a validated read or fragment"
         strand(rec) == featurestrand || continue
diff --git a/gff_feature.jl b/src/gff_feature.jl
similarity index 100%
rename from gff_feature.jl
rename to src/gff_feature.jl
diff --git a/test_data/tiny/expected.counts.txt b/test/data/expected.counts.txt
similarity index 100%
rename from test_data/tiny/expected.counts.txt
rename to test/data/expected.counts.txt
diff --git a/test_data/tiny/tiny.bam b/test/data/tiny.bam
similarity index 100%
rename from test_data/tiny/tiny.bam
rename to test/data/tiny.bam
diff --git a/test_data/tiny/tiny.bam.bai b/test/data/tiny.bam.bai
similarity index 100%
rename from test_data/tiny/tiny.bam.bai
rename to test/data/tiny.bam.bai
diff --git a/test/data/tiny.exclude.bed b/test/data/tiny.exclude.bed
new file mode 100644
index 0000000..f572301
--- /dev/null
+++ b/test/data/tiny.exclude.bed
@@ -0,0 +1 @@
+NZ_CP059840.2	4595849	4596275	exclude_cds	0	+
diff --git a/test_data/tiny/tiny.gff3 b/test/data/tiny.gff3
similarity index 100%
rename from test_data/tiny/tiny.gff3
rename to test/data/tiny.gff3
diff --git a/test_data/tiny/run_test.sh b/test/run_test.sh
similarity index 59%
rename from test_data/tiny/run_test.sh
rename to test/run_test.sh
index c25a021..2632cbc 100755
--- a/test_data/tiny/run_test.sh
+++ b/test/run_test.sh
@@ -3,16 +3,16 @@
 # Override the FADU command with FADU_CMD (e.g. FADU_CMD="fadu.jl" for an installed copy).
 set -euo pipefail
 
-TEST_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
-ROOT_DIR="$(cd "${TEST_DIR}/../.." && pwd)"
-FADU_CMD="${FADU_CMD:-julia --project=${ROOT_DIR}/fadu_pkgs ${ROOT_DIR}/fadu.jl}"
+DATA_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)/data"
+ROOT_DIR="$(cd "${DATA_DIR}/../.." && pwd)"
+FADU_CMD="${FADU_CMD:-julia ${ROOT_DIR}/fadu.jl}"
 
 OUT_DIR="$(mktemp -d)"
 trap 'rm -rf "${OUT_DIR}"' EXIT
 
-${FADU_CMD} -g "${TEST_DIR}/tiny.gff3" -b "${TEST_DIR}/tiny.bam" -o "${OUT_DIR}" -s "reverse" -f "CDS" -a "ID"
+${FADU_CMD} -g "${DATA_DIR}/tiny.gff3" -b "${DATA_DIR}/tiny.bam" -o "${OUT_DIR}" -s "reverse" -f "CDS" -a "ID"
 
-if diff "${TEST_DIR}/expected.counts.txt" "${OUT_DIR}/tiny.counts.txt"; then
+if diff "${DATA_DIR}/expected.counts.txt" "${OUT_DIR}/tiny.counts.txt"; then
     echo "PASS: output matches expected.counts.txt"
 else
     echo "FAIL: output differs from expected.counts.txt" >&2
diff --git a/test/runtests.jl b/test/runtests.jl
new file mode 100644
index 0000000..6912efe
--- /dev/null
+++ b/test/runtests.jl
@@ -0,0 +1,67 @@
+using FADU
+using Logging
+using Test
+
+const DATA_DIR = joinpath(@__DIR__, "data")
+const TINY_ARGS = ["-g", joinpath(DATA_DIR, "tiny.gff3"), "-b", joinpath(DATA_DIR, "tiny.bam"),
+                   "-f", "CDS", "-a", "ID"]
+const EXCLUDE_BED = joinpath(DATA_DIR, "tiny.exclude.bed")
+const EXCLUDED_FEATURE = "cds-WP_001339389.1"    # Entirely covered by tiny.exclude.bed
+
+# Run FADU quietly and return the exit code and the lines of the counts file
+function run_fadu(extra_args::Vector{String}=String[]; stranded::String="reverse")
+    mktempdir() do out_dir
+        exit_code = with_logger(NullLogger()) do
+            FADU.main([TINY_ARGS; "-s"; stranded; "-o"; out_dir; extra_args])
+        end
+        return exit_code, readlines(joinpath(out_dir, "tiny.counts.txt"))
+    end
+end
+
+@testset "FADU" begin
+    @testset "version matches Project.toml" begin
+        @test FADU.VERSION_NUMBER == string(pkgversion(FADU))
+    end
+
+    expected = readlines(joinpath(DATA_DIR, "expected.counts.txt"))
+
+    @testset "tiny dataset matches expected output" begin
+        exit_code, lines = run_fadu()
+        @test exit_code == 0
+        @test lines == expected
+    end
+
+    @testset "--no_output_header omits header line" begin
+        exit_code, lines = run_fadu(["--no_output_header"])
+        @test exit_code == 0
+        @test lines == expected[2:end]
+    end
+
+    @testset "--remove_multimapped runs" begin
+        exit_code, lines = run_fadu(["--remove_multimapped"])
+        @test exit_code == 0
+        @test length(lines) == length(expected)
+    end
+
+    @testset "--exclude_regions_file (stranded=$stranded)" for stranded in ["no", "reverse"]
+        _, baseline = run_fadu(; stranded)
+        exit_code, lines = run_fadu(["--exclude_regions_file", EXCLUDE_BED]; stranded)
+        @test exit_code == 0
+        # The excluded feature loses all of its unique coordinates and counts...
+        excluded = split(only(filter(l -> startswith(l, EXCLUDED_FEATURE * "\t"), lines)), '\t')
+        @test excluded[2] == "0"
+        @test excluded[4] == "0.00"
+        # ...while every other feature keeps its unique coordinates
+        uniq_lens(ls) = Dict(f[1] => f[2] for f in split.(ls[2:end], '\t') if f[1] != EXCLUDED_FEATURE)
+        @test uniq_lens(lines) == uniq_lens(baseline)
+    end
+
+    @testset "missing BAM file returns non-zero exit code" begin
+        mktempdir() do out_dir
+            args = ["-g", joinpath(DATA_DIR, "tiny.gff3"), "-b", joinpath(DATA_DIR, "missing.bam"), "-o", out_dir]
+            @test with_logger(NullLogger()) do
+                FADU.main(args)
+            end == 1
+        end
+    end
+end