From 1bbdeb15ef36bff0292636a9978ff28c8a2139b4 Mon Sep 17 00:00:00 2001 From: adkinsrs Date: Wed, 28 Oct 2026 08:56:53 -0400 Subject: [PATCH] Fix BAM index lookup and --log_level handling (v1.10.1) - Only replace the final ".bam" extension when looking for .bai, then fall back to .bam.bai. Previously a BAM path without a ".bam" extension was used as its own index, and ".bam" in directory names was rewritten. A missing index is now reported during argument validation. - Make --log_level case-insensitive and drop CRITICAL (never accepted) from the help text. - Show "fadu" as the program name and fix the --feature_type and --attribute_type help text; regenerate the usage block in the docs. - Add tests for the index lookup, extensionless BAMs, a missing index, and lowercase log levels. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 6 ++++++ Manifest.toml | 2 +- Project.toml | 2 +- README.md | 33 +++++++++++++++++---------------- docs/index.md | 33 +++++++++++++++++---------------- src/FADU.jl | 41 ++++++++++++++++++++--------------------- test/runtests.jl | 48 ++++++++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 110 insertions(+), 55 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3bae9c8..4ee43f8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,11 @@ # CHANGELOG +## v1.10.1 + +* BUGFIX - BAM index lookup replaced every ".bam" in the path with ".bai". A BAM path without a ".bam" extension was treated as its own index (failing with "input is not a valid BAI file"), and ".bam" in a directory name was rewritten. Now only the final ".bam" extension is replaced, then `.bam.bai` is tried. A missing index is reported during argument validation. +* BUGFIX - `--log_level` is now case-insensitive, and CRITICAL (which was never accepted) is removed from the help text. +* Help text now shows `fadu` as the program name and fixes the `--feature_type`/`--attribute_type` descriptions (GTF is not supported; column 9 is an attribute). + ## 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. diff --git a/Manifest.toml b/Manifest.toml index 2d10f1e..4b25589 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -127,7 +127,7 @@ version = "1.11.0" deps = ["ArgParse", "BED", "BGZFStreams", "GFF3", "GenomicFeatures", "Indexes", "Logging", "PrecompileTools", "Printf", "StructArrays", "XAM"] path = "." uuid = "866296ad-388f-4871-8691-d7485d094d7e" -version = "1.10.0" +version = "1.10.1" [[deps.FASTX]] deps = ["Automa", "BioGenerics", "PrecompileTools", "ScanByte", "StringViews", "TranscodingStreams"] diff --git a/Project.toml b/Project.toml index a7e61cf..057534c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FADU" uuid = "866296ad-388f-4871-8691-d7485d094d7e" -version = "1.10.0" +version = "1.10.1" authors = ["adkinsrs "] [deps] diff --git a/README.md b/README.md index d8143d3..0b3e3eb 100644 --- a/README.md +++ b/README.md @@ -17,7 +17,7 @@ Most current available quantification tools for transcriptomics analyses have be ## Current FADU release -v1.10.0 +v1.10.1 ## Requirements @@ -137,12 +137,12 @@ By default the script runs `julia ./fadu.jl`. Set the `FADU_CMD` environment va ## Usage ```bash -usage: fadu.jl -b /path/to/file.bam -g /path/to/annotation.gff3 - -o /path/to/output/dir/ [--no_output_header] - [-x /path/to/regions.bed] [-s STRANDED] - [-f FEATURE_TYPE] [-a ATTRIBUTE_TYPE] [-p] - [-m MAX_FRAGMENT_SIZE] [-M] [-e EM_ITER] [--version] - [-h] +usage: fadu -b /path/to/file.bam -g /path/to/annotation.gff3 + -o /path/to/output/dir/ [--no_output_header] + [-x /path/to/regions.bed] [-s STRANDED] [-f FEATURE_TYPE] + [-a ATTRIBUTE_TYPE] [-p] [-m MAX_FRAGMENT_SIZE] [-M] + [-e EM_ITER] [--log_level LEVEL] + [--log_file /path/to/logfile.log] [--version] [-h] Generate counts of reads that map to non-overlapping portions of genes @@ -165,13 +165,13 @@ optional arguments: strand-specific assay. Choose between 'yes', 'no', or 'reverse'. (default: "no") -f, --feature_type FEATURE_TYPE - Which GFF3/GTF feature type (column 3) to - obtain. readcount statistics for. - Case-sensitive. (default: "gene") + Which GFF3 feature type (column 3) to obtain + readcount statistics for. Case-sensitive. + (default: "gene") -a, --attribute_type ATTRIBUTE_TYPE - Which GFF3/GTF feature type (column 9) to - obtain. readcount statistics for. - Case-sensitive. (default: "ID") + Which GFF3 attribute (column 9) to use as the + feature ID in the output. Case-sensitive. + (default: "ID") -p, --keep_only_proper_pairs If enabled, keep only properly paired reads when performing calculations. @@ -187,7 +187,7 @@ optional arguments: -M, --remove_multimapped If enabled, remove any reads or fragments that are mapped to multiple regions of the genome, - indiated by the 'NH' attribute being greater + indicated by the 'NH' attribute being greater than 1. Otherwise, use EM algorithm to quantify reads after all other reads are counted. @@ -196,8 +196,9 @@ optional arguments: on multimapped read depth. Only applies if --remove_multimapped flag is not passed (is disabled) (type: Int64, default: 1) - --log_level LEVEL Set the log level. Options are: DEBUG, INFO, - WARNING, ERROR, CRITICAL. (default: "INFO") + --log_level LEVEL Set the log level (case-insensitive). Options + are: DEBUG, INFO, WARNING, ERROR. (default: + "INFO") --log_file /path/to/logfile.log Path to log file. If not specified, log messages will be printed to STDERR. diff --git a/docs/index.md b/docs/index.md index 5e6450a..38fc453 100644 --- a/docs/index.md +++ b/docs/index.md @@ -17,7 +17,7 @@ Most current available quantification tools for transcriptomics analyses have be ## Current FADU release -v1.10.0 +v1.10.1 ## Requirements @@ -138,12 +138,12 @@ By default the script runs `julia ./fadu.jl`. Set the `FADU_CMD` environment va ## Usage ```bash -usage: fadu.jl -b /path/to/file.bam -g /path/to/annotation.gff3 - -o /path/to/output/dir/ [--no_output_header] - [-x /path/to/regions.bed] [-s STRANDED] - [-f FEATURE_TYPE] [-a ATTRIBUTE_TYPE] [-p] - [-m MAX_FRAGMENT_SIZE] [-M] [-e EM_ITER] [--version] - [-h] +usage: fadu -b /path/to/file.bam -g /path/to/annotation.gff3 + -o /path/to/output/dir/ [--no_output_header] + [-x /path/to/regions.bed] [-s STRANDED] [-f FEATURE_TYPE] + [-a ATTRIBUTE_TYPE] [-p] [-m MAX_FRAGMENT_SIZE] [-M] + [-e EM_ITER] [--log_level LEVEL] + [--log_file /path/to/logfile.log] [--version] [-h] Generate counts of reads that map to non-overlapping portions of genes @@ -166,13 +166,13 @@ optional arguments: strand-specific assay. Choose between 'yes', 'no', or 'reverse'. (default: "no") -f, --feature_type FEATURE_TYPE - Which GFF3/GTF feature type (column 3) to - obtain. readcount statistics for. - Case-sensitive. (default: "gene") + Which GFF3 feature type (column 3) to obtain + readcount statistics for. Case-sensitive. + (default: "gene") -a, --attribute_type ATTRIBUTE_TYPE - Which GFF3/GTF feature type (column 9) to - obtain. readcount statistics for. - Case-sensitive. (default: "ID") + Which GFF3 attribute (column 9) to use as the + feature ID in the output. Case-sensitive. + (default: "ID") -p, --keep_only_proper_pairs If enabled, keep only properly paired reads when performing calculations. @@ -188,7 +188,7 @@ optional arguments: -M, --remove_multimapped If enabled, remove any reads or fragments that are mapped to multiple regions of the genome, - indiated by the 'NH' attribute being greater + indicated by the 'NH' attribute being greater than 1. Otherwise, use EM algorithm to quantify reads after all other reads are counted. @@ -197,8 +197,9 @@ optional arguments: on multimapped read depth. Only applies if --remove_multimapped flag is not passed (is disabled) (type: Int64, default: 1) - --log_level LEVEL Set the log level. Options are: DEBUG, INFO, - WARNING, ERROR, CRITICAL. (default: "INFO") + --log_level LEVEL Set the log level (case-insensitive). Options + are: DEBUG, INFO, WARNING, ERROR. (default: + "INFO") --log_file /path/to/logfile.log Path to log file. If not specified, log messages will be printed to STDERR. diff --git a/src/FADU.jl b/src/FADU.jl index e20f0d7..8da2a45 100644 --- a/src/FADU.jl +++ b/src/FADU.jl @@ -40,7 +40,7 @@ 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 VERSION_NUMBER = "1.10.1" # 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 @@ -100,7 +100,7 @@ end 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", + prog = "fadu", version = VERSION_NUMBER, add_version = true, add_help = true) @@ -131,10 +131,10 @@ function parse_commandline(cli_args::AbstractVector{<:AbstractString}) 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." + help = "Which GFF3 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." + help = "Which GFF3 attribute (column 9) to use as the feature ID in the output. Case-sensitive." default = "ID" "--keep_only_proper_pairs", "-p" help = "If enabled, keep only properly paired reads when performing calculations." @@ -146,7 +146,7 @@ function parse_commandline(cli_args::AbstractVector{<:AbstractString}) 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." + help = "If enabled, remove any reads or fragments that are mapped to multiple regions of the genome, indicated 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" @@ -156,15 +156,13 @@ function parse_commandline(cli_args::AbstractVector{<:AbstractString}) range_tester = (x->x>0) dest_name = "em_iter" "--log_level" - help = "Set the log level. Options are: DEBUG, INFO, WARNING, ERROR, CRITICAL." + help = "Set the log level (case-insensitive). Options are: DEBUG, INFO, WARNING, ERROR." default = "INFO" metavar = "LEVEL" - range_tester = (x->x in ["DEBUG", "INFO", "WARNING", "ERROR"]) + range_tester = (x->uppercase(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) @@ -181,11 +179,23 @@ function setup_logger(args) end +function find_bam_index(bam_file::AbstractString) + """Return the path to the index of the BAM file (.bai or .bam.bai), or nothing if neither exists.""" + root, ext = splitext(bam_file) + candidates = lowercase(ext) == ".bam" ? [root * ".bai", bam_file * ".bai"] : [bam_file * ".bai"] + for candidate in candidates + isfile(candidate) && return candidate + end + return nothing +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") + args["bam_index"] = find_bam_index(args["bam_file"]) + args["bam_index"] === nothing && error("Could not find an index file for BAM file $(args["bam_file"]). Expected .bai or .bam.bai in the same directory") if !isdir(args["output_dir"]) @info("Creating output directory at ", args["output_dir"]) mkdir(args["output_dir"]) @@ -242,20 +252,9 @@ function main(cli_args::AbstractVector{<:AbstractString}=ARGS) @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) + reader = open(BAM.Reader, args["bam_file"], index = args["bam_index"]) @info("Now finding overlaps between alignment and annotation records...") multimapped_dict = process_all_feature_overlaps(feat_overlaps, features, reader, args) diff --git a/test/runtests.jl b/test/runtests.jl index 6912efe..2e97aff 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -56,6 +56,54 @@ end @test uniq_lens(lines) == uniq_lens(baseline) end + @testset "BAM index lookup" begin + mktempdir() do dir + # Extensionless BAM must not be mistaken for its own index + bam = joinpath(dir, "sample") + cp(joinpath(DATA_DIR, "tiny.bam"), bam) + @test FADU.find_bam_index(bam) === nothing + cp(joinpath(DATA_DIR, "tiny.bam.bai"), bam * ".bai") + @test FADU.find_bam_index(bam) == bam * ".bai" + + # Only the final ".bam" extension is replaced + sub = mkpath(joinpath(dir, "run.bam_files")) + bam = joinpath(sub, "x.bam") + cp(joinpath(DATA_DIR, "tiny.bam"), bam) + cp(joinpath(DATA_DIR, "tiny.bam.bai"), joinpath(sub, "x.bai")) + @test FADU.find_bam_index(bam) == joinpath(sub, "x.bai") + end + end + + @testset "extensionless BAM runs" begin + mktempdir() do dir + bam = joinpath(dir, "sample") + cp(joinpath(DATA_DIR, "tiny.bam"), bam) + cp(joinpath(DATA_DIR, "tiny.bam.bai"), bam * ".bai") + args = ["-g", joinpath(DATA_DIR, "tiny.gff3"), "-b", bam, "-o", dir, "-s", "reverse", "-f", "CDS", "-a", "ID"] + @test with_logger(NullLogger()) do + FADU.main(args) + end == 0 + @test readlines(joinpath(dir, "sample.counts.txt")) == expected + end + end + + @testset "missing BAM index returns non-zero exit code" begin + mktempdir() do dir + bam = joinpath(dir, "tiny.bam") + cp(joinpath(DATA_DIR, "tiny.bam"), bam) + args = ["-g", joinpath(DATA_DIR, "tiny.gff3"), "-b", bam, "-o", dir] + @test with_logger(NullLogger()) do + FADU.main(args) + end == 1 + end + end + + @testset "--log_level is case-insensitive" begin + exit_code, lines = run_fadu(["--log_level", "warning"]) + @test exit_code == 0 + @test lines == expected + 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]