Fix BAM index lookup and --log_level handling (v1.10.1)

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 `<name>.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.

Project.toml

@@ -1,6 +1,6 @@
 name = "FADU"
 uuid = "866296ad-388f-4871-8691-d7485d094d7e"
-version = "1.10.0"
+version = "1.10.1"
 authors = ["adkinsrs <sadkins@som.umaryland.edu>"]
 
 [deps]

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.

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.

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 (<name>.bai or <name>.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 <name>.bai or <name>.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)

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]