From 0ed2cc4090ebad8fe418f9324529de40f3c8a3a1 Mon Sep 17 00:00:00 2001 From: DanieCuevas <43822444+DanielCuevas1208@users.noreply.github.com> Date: Sun, 2 Aug 2026 19:16:08 -0700 Subject: [PATCH] feat: extend codon atlas --- .github/workflows/ci.yml | 11 +- CHANGELOG.md | 18 +++ README.md | 126 ++++++++++++-------- codonatlas.nimble | 4 +- examples/demo_fastq.nim | 11 ++ src/codonatlas.nim | 9 +- src/codonatlas/fastq.nim | 97 ++++++++++++++++ src/codonatlas/quality.nim | 50 ++++++++ src/codonatlas/report.nim | 99 +++++++++++++--- src/codonatlas_cli.nim | 59 +++++++--- tests/data/expected_fastq_report.txt | 39 +++++++ tests/data/sample.fastq | 8 ++ tests/test_fastq.nim | 166 +++++++++++++++++++++++++++ 13 files changed, 614 insertions(+), 83 deletions(-) create mode 100644 examples/demo_fastq.nim create mode 100644 src/codonatlas/fastq.nim create mode 100644 src/codonatlas/quality.nim create mode 100644 tests/data/expected_fastq_report.txt create mode 100644 tests/data/sample.fastq create mode 100644 tests/test_fastq.nim diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 884a626..dffb5bb 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -33,7 +33,7 @@ jobs: - name: Build release CLI run: nim c -d:release --path:src -o:build/codonatlas src/codonatlas_cli.nim - - name: Smoke test the CLI + - name: Smoke test the FASTA path shell: bash run: | BIN="build/codonatlas" @@ -41,3 +41,12 @@ jobs: BIN="${BIN}.exe" fi "./${BIN}" tests/data/sample.fasta >/dev/null + + - name: Smoke test the FASTQ path + shell: bash + run: | + BIN="build/codonatlas" + if [[ "$RUNNER_OS" == "Windows" ]]; then + BIN="${BIN}.exe" + fi + "./${BIN}" tests/data/sample.fastq --quality=20 --kmer=3 >/dev/null diff --git a/CHANGELOG.md b/CHANGELOG.md index d845630..8a50c4e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,24 @@ This project follows semantic versioning. +## [0.2.0] - 2026-08-03 + +Second coherent release. This release completes the FASTQ milestone from the +roadmap. + +Added: + +- FASTQ parsing with headers, descriptions, wrapped lines, and quality lines. +- Phred+33 quality score statistics per record. +- A FASTQ text report with per-record quality summaries. +- A `--format` option and automatic format detection by file extension. +- A `--quality` option to add a below-threshold count line. +- A bundled FASTQ fixture and a golden report comparison. +- Unit tests for the FASTQ parser and the quality statistics. + +The package uses only the Nim standard library. It has no external +dependencies, so it needs no lockfile. + ## [0.1.0] - 2026-08-02 First coherent release. diff --git a/README.md b/README.md index 05fcee3..1b2b562 100644 --- a/README.md +++ b/README.md @@ -1,15 +1,18 @@ # Codon Atlas -Codon Atlas is a Nim library for DNA and RNA sequence files. It computes -composition statistics. These include GC content, k-mer counts, and codon -usage tables. A command-line tool writes the results to a text report. +Codon Atlas is a Nim library for DNA and RNA sequence files. It reads FASTA +and FASTQ files. It computes GC content, k-mer counts, and codon usage tables. +It summarizes FASTQ quality scores. A command-line tool writes the results to +a text report. ## Features - Read FASTA files with a permissive parser. +- Read FASTQ files with their quality lines. - Compute GC content and per-base counts. - Count overlapping k-mers. - Build codon usage tables from the standard genetic code. +- Summarize FASTQ quality scores per record. - Write a text report to the terminal or a file. - Share one API across the library and the CLI. @@ -31,6 +34,7 @@ Run the test suite: ``` nim c -r --path:src --outdir:build/tests tests/test_codonatlas.nim +nim c -r --path:src --outdir:build/tests tests/test_fastq.nim ``` On a checked-out repository, `nimble build` and `nimble test` are the @@ -64,6 +68,25 @@ Show only the top 10 k-mers: ./build/codonatlas tests/data/sample.fasta --kmer=3 --top=10 ``` +Produce a report for a FASTQ file. The tool detects the format from the file +extension: + +``` +./build/codonatlas tests/data/sample.fastq +``` + +Add a quality threshold line to a FASTQ report: + +``` +./build/codonatlas tests/data/sample.fastq --quality=20 +``` + +Force a format when the file extension is not standard: + +``` +./build/codonatlas reads.data --format=fastq +``` + ### Library Import the package and build a report in code: @@ -75,70 +98,69 @@ let records = readFasta("tests/data/sample.fasta") echo buildReport(records, source = "sample.fasta") ``` -Run the bundled example: +Build a report from a FASTQ file: + +```nim +import codonatlas + +let records = readFastq("tests/data/sample.fastq") +echo buildFastqReport(records, source = "sample.fastq", + options = ReportOptions(qualityThreshold = 20)) +``` + +Run the bundled examples: ``` nim c -r --path:src --outdir:build/examples examples/demo.nim +nim c -r --path:src --outdir:build/examples examples/demo_fastq.nim ``` ## Sample output -The command below prints a report for the bundled sample file: +The command below prints a FASTQ report for the bundled sample file: ``` -./build/codonatlas tests/data/sample.fasta --kmer=3 --top=3 +./build/codonatlas tests/data/sample.fastq --quality=20 ``` ``` Codon Atlas report ================== -Source: tests/data/sample.fasta -Generated: 2026-08-03T00:27:23Z +Source: tests/data/sample.fastq +Generated: 2026-08-03T02:15:15Z Records: 2 -Total canonical bases: 181 +Total canonical bases: 35 Nucleotide composition ---------------------- -Record 1: demo_orf - Description: synthetic coding sequence, 144 bp, starts with ATG and ends with TAA - Length: 144 canonical bases - Bases: A 40 C 31 G 38 T 35 U 0 N 0 - GC content: 47.9% - -Record 2: promoter_fragment - Description: synthetic non-coding fragment - Length: 37 canonical bases - Bases: A 9 C 8 G 9 T 11 U 0 N 0 - GC content: 45.9% +Record 1: read1 + Description: synthetic high-quality read, 23 bp + Length: 23 canonical bases + Bases: A 5 C 5 G 7 T 6 U 0 N 0 + GC content: 52.2% + Quality: mean 38.3, min 37, max 40 + Bases below Q20: 0 + +Record 2: read2 + Description: synthetic low-quality read, 12 bp + Length: 12 canonical bases + Bases: A 3 C 4 G 2 T 3 U 0 N 0 + GC content: 50.0% + Quality: mean 7.0, min 0, max 13 + Bases below Q20: 12 Overall ------- -Total canonical bases: 181 -GC content: 47.5% - -K-mer counts (k = 3) --------------------------------- -Top 3 of 61 distinct k-mers. -CGT 8 -GTA 7 -GCG 6 - -Codon usage table (reading frame 0) ------------------------------------ -Total codons: 60 - -Count Frequency Per 1000 Codon AA Amino acid --------------------------------------------- - 4 0.0667 66.7 CGU R Arginine - 3 0.0500 50.0 AUG M Methionine - 3 0.0500 50.0 GAC D Aspartic acid - 3 0.0500 50.0 GAU D Aspartic acid - 2 0.0333 33.3 AAG K Lysine +Total canonical bases: 35 +GC content: 51.4% ``` -The report lists all 64 codons of the standard genetic code. The output above -is abbreviated. Codons use RNA letters (U instead of T). A stop codon shows as -`*`. +Quality scores use the Phred+33 encoding. A higher score means a more +confident base call. A score of 20 is a common quality threshold. + +The FASTA report adds a codon usage table. It lists all 64 codons of the +standard genetic code. Codons use RNA letters (U instead of T). A stop codon +shows as `*`. ## Project layout @@ -146,12 +168,15 @@ is abbreviated. Codons use RNA letters (U instead of T). A stop codon shows as codonatlas.nimble package metadata and tasks src/codonatlas.nim public API entry point src/codonatlas/fasta.nim FASTA parsing +src/codonatlas/fastq.nim FASTQ parsing +src/codonatlas/quality.nim quality score statistics src/codonatlas/composition.nim GC content, base counts, k-mers src/codonatlas/gencode.nim standard genetic code src/codonatlas/codons.nim codon usage tables -src/codonatlas/report.nim text report builder +src/codonatlas/report.nim text report builders src/codonatlas_cli.nim command-line interface -examples/demo.nim runnable example +examples/demo.nim runnable FASTA example +examples/demo_fastq.nim runnable FASTQ example tests/ unit tests and fixtures ``` @@ -159,7 +184,7 @@ tests/ unit tests and fixtures Run `nimble check` to typecheck every module. Run `nimble test` to execute the unit tests. The suite covers parsing, statistics, the genetic code, the -report builder, and a golden report comparison. +report builders, quality scores, and golden report comparisons. ## Limitations @@ -167,10 +192,15 @@ report builder, and a golden report comparison. - Codon counts use reading frame 0 on the given strands. - Ambiguous bases (N and others) are excluded from composition ratios. - The table normalizes codons to RNA form. It merges DNA and RNA input. +- The FASTQ parser assumes Sanger quality scores. It rejects sequences and + quality lines with different lengths. +- FASTQ reports do not include codon usage tables. Reads are usually short + unassembled fragments. ## Roadmap -Release 2 adds FASTQ support. It will summarize quality scores per record. +Release 2 is complete. Codon Atlas reads FASTQ files and summarizes quality +scores per record. Release 3 adds translation. It will map codons to amino acids. It will also detect open reading frames. diff --git a/codonatlas.nimble b/codonatlas.nimble index 5d5e7c4..70461e5 100644 --- a/codonatlas.nimble +++ b/codonatlas.nimble @@ -1,4 +1,4 @@ -version = "0.1.0" +version = "0.2.0" author = "Codon Atlas contributors" description = "Composition statistics for DNA and RNA sequence files" license = "MIT" @@ -10,8 +10,10 @@ requires "nim >= 1.6.0" task test, "Run the unit test suite": exec "nim c -r --path:src --outdir:build/tests tests/test_codonatlas.nim" + exec "nim c -r --path:src --outdir:build/tests tests/test_fastq.nim" task check, "Typecheck all modules": exec "nim check --path:src src/codonatlas.nim" exec "nim check --path:src src/codonatlas_cli.nim" exec "nim check --path:src tests/test_codonatlas.nim" + exec "nim check --path:src tests/test_fastq.nim" diff --git a/examples/demo_fastq.nim b/examples/demo_fastq.nim new file mode 100644 index 0000000..5432d7b --- /dev/null +++ b/examples/demo_fastq.nim @@ -0,0 +1,11 @@ +## Runnable example: build a Codon Atlas report for the bundled FASTQ file. + +import std/os +import codonatlas + +let samplePath = currentSourcePath().parentDir() / ".." / "tests" / "data" / "sample.fastq" +let records = readFastq(samplePath) + +stdout.write buildFastqReport(records, source = "sample.fastq", + options = ReportOptions(kmer: 3, qualityThreshold: 20), + generatedAt = "2026-08-02T12:00:00Z") diff --git a/src/codonatlas.nim b/src/codonatlas.nim index 77e0c45..f2fcb0d 100644 --- a/src/codonatlas.nim +++ b/src/codonatlas.nim @@ -1,8 +1,9 @@ ## Codon Atlas: composition statistics for DNA and RNA sequence files. ## -## The package reads FASTA files, computes GC content and k-mer counts, and -## builds codon usage tables. A CLI writes the results to a plain-text report. +## The package reads FASTA and FASTQ files, computes GC content and k-mer +## counts, and builds codon usage tables. It summarizes FASTQ quality scores. +## A CLI writes the results to a plain-text report. -import codonatlas/[fasta, composition, gencode, codons, report] +import codonatlas/[fasta, fastq, quality, composition, gencode, codons, report] -export fasta, composition, gencode, codons, report +export fasta, fastq, quality, composition, gencode, codons, report diff --git a/src/codonatlas/fastq.nim b/src/codonatlas/fastq.nim new file mode 100644 index 0000000..3a925a7 --- /dev/null +++ b/src/codonatlas/fastq.nim @@ -0,0 +1,97 @@ +## FASTQ reading and parsing. +## +## The parser accepts the standard Sanger FASTQ format. Each record has a +## header line, a sequence, a `+` separator, and a quality line. Sequence and +## quality lines may wrap over multiple lines. Their combined lengths must +## match. Quality scores use the Phred+33 encoding. + +import std/strutils + +type + FastqError* = object of CatchableError + ## Raised when the FASTQ text is malformed. + + FastqRecord* = object + id*: string + ## First token of the header line, after the `@`. + description*: string + ## Remainder of the header line, if any. + sequence*: string + ## Concatenated sequence bases, with whitespace removed. + quality*: string + ## Concatenated quality characters, with whitespace removed. + +proc stripWhitespace(line: string): string = + result = newStringOfCap(line.len) + for ch in line: + if ch notin {' ', '\t', '\r', '\n'}: + result.add ch + +proc parseHeader(line: string, lineNo: int): tuple[id, description: string] = + let header = line[1 .. ^1].strip + if header.len == 0: + raise newException(FastqError, "empty sequence header on line " & $lineNo) + let fields = header.splitWhitespace() + result.id = fields[0] + result.description = fields[1 .. ^1].join(" ") + +proc parseFastq*(content: string): seq[FastqRecord] = + ## Parse FASTQ text into a list of records. + ## + ## Raises `FastqError` when the text is malformed. An empty input produces + ## an empty list. + let lines = splitLines(content) + var index = 0 + while index < lines.len: + if lines[index].strip.len == 0: + inc index + continue + if lines[index][0] != '@': + raise newException(FastqError, + "expected a header starting with '@' (line " & $(index + 1) & ")") + var record = FastqRecord() + let (id, description) = parseHeader(lines[index], index + 1) + record.id = id + record.description = description + inc index + + var sequenceParts: seq[string] + while index < lines.len: + let line = lines[index] + if line.len > 0 and line[0] == '+': + break + if line.len > 0: + sequenceParts.add stripWhitespace(line) + inc index + if sequenceParts.len == 0: + raise newException(FastqError, + "record '" & record.id & "' has no sequence") + if index >= lines.len: + raise newException(FastqError, + "record '" & record.id & "' is missing the '+' separator") + record.sequence = sequenceParts.join("") + inc index + + if index >= lines.len: + raise newException(FastqError, + "record '" & record.id & "' has no quality line") + var quality = "" + while index < lines.len and quality.len < record.sequence.len: + if lines[index].len > 0: + quality.add stripWhitespace(lines[index]) + inc index + if quality.len < record.sequence.len: + raise newException(FastqError, + "quality is shorter than the sequence in record '" & record.id & "'") + if quality.len > record.sequence.len: + raise newException(FastqError, + "quality is longer than the sequence in record '" & record.id & "'") + record.quality = quality + result.add record + +proc readFastq*(path: string): seq[FastqRecord] = + ## Read a FASTQ file and parse its records. + ## + ## Raises `IOError` when the file cannot be read. Raises `FastqError` when + ## the file content is malformed. + parseFastq(readFile(path)) diff --git a/src/codonatlas/quality.nim b/src/codonatlas/quality.nim new file mode 100644 index 0000000..062a0c2 --- /dev/null +++ b/src/codonatlas/quality.nim @@ -0,0 +1,50 @@ +## Quality score statistics for FASTQ sequences. +## +## Scores use the Sanger Phred+33 encoding. Each quality character maps to a +## score of `ord(c) - 33`. A valid score is an integer from 0 to 93. + +const StandardQualityOffset* = 33 + ## Phred offset for Sanger FASTQ quality strings. + +type + QualitySummary* = object + count*: int + ## Number of quality scores. + min*: int + ## Lowest score. + max*: int + ## Highest score. + mean*: float + ## Average score. Zero when there are no scores. + +proc qualityScore*(base: char, offset = StandardQualityOffset): int = + ## Return the Phred quality score for a base character. + ## + ## Raises `ValueError` when the score falls outside 0 to 93. + result = ord(base) - offset + if result < 0 or result > 93: + raise newException(ValueError, "quality score out of range: " & $result) + +proc qualitySummary*(quality: string, offset = StandardQualityOffset): QualitySummary = + ## Summarize a quality string. + ## + ## An empty string produces a summary with all values at zero. + if quality.len == 0: + return + result.count = quality.len + result.min = high(int) + result.max = low(int) + var total = 0 + for base in quality: + let score = qualityScore(base, offset) + total += score + result.min = min(result.min, score) + result.max = max(result.max, score) + result.mean = total.float / result.count.float + +proc basesBelowQuality*(quality: string, threshold: int, + offset = StandardQualityOffset): int = + ## Count positions with a score strictly below `threshold`. + for base in quality: + if qualityScore(base, offset) < threshold: + inc result diff --git a/src/codonatlas/report.nim b/src/codonatlas/report.nim index 928e86f..d207af7 100644 --- a/src/codonatlas/report.nim +++ b/src/codonatlas/report.nim @@ -9,6 +9,8 @@ import std/times import codonatlas/codons import codonatlas/composition import codonatlas/fasta +import codonatlas/fastq +import codonatlas/quality type ReportOptions* = object @@ -16,6 +18,9 @@ type ## When greater than 0, add a k-mer count section for this k. kmerTop*: int ## When greater than 0, show only the top k-mers in that section. + qualityThreshold*: int + ## When greater than 0, add a below-threshold count for each FASTQ + ## record's quality line. proc cmpKmer(a, b: tuple[key: string, val: int]): int = if a.val != b.val: @@ -30,6 +35,22 @@ proc formatCounts(counts: NucleotideCounts): string = "A " & $counts.a & " C " & $counts.c & " G " & $counts.g & " T " & $counts.t & " U " & $counts.u & " N " & $counts.n +proc kmerSection(sequences: openArray[string], k, top: int): seq[string] = + var combined = initCountTable[string]() + for sequence in sequences: + for key, value in kmerCounts(sequence, k): + combined.inc(key, value) + var items = combined.pairs.toSeq + items.sort(cmpKmer) + let shown = if top > 0: min(top, items.len) else: items.len + result.add "" + result.add "K-mer counts (k = " & $k & ")" + result.add "--------------------------------" + if top > 0: + result.add "Top " & $shown & " of " & $items.len & " distinct k-mers." + for i in 0 ..< shown: + result.add alignLeft(items[i][0], 4) & " " & align($items[i][1], 6) + proc gcPercentOf(counts: NucleotideCounts): float = let total = totalBases(counts) if total == 0: @@ -86,20 +107,7 @@ proc buildReport*(records: seq[Record], source: string, outLines.add "GC content: " & fmtDecimal(gcPercentOf(allCounts), 1) & "%" if options.kmer > 0: - var combined = initCountTable[string]() - for record in records: - for key, value in kmerCounts(record.sequence, options.kmer): - combined.inc(key, value) - var items = combined.pairs.toSeq - items.sort(cmpKmer) - let shown = if options.kmerTop > 0: min(options.kmerTop, items.len) else: items.len - outLines.add "" - outLines.add "K-mer counts (k = " & $options.kmer & ")" - outLines.add "--------------------------------" - if options.kmerTop > 0: - outLines.add "Top " & $shown & " of " & $items.len & " distinct k-mers." - for i in 0 ..< shown: - outLines.add alignLeft(items[i][0], 4) & " " & align($items[i][1], 6) + outLines.add kmerSection(records.mapIt(it.sequence), options.kmer, options.kmerTop) let usage = countCodonUsageMany(records.mapIt(it.sequence)) outLines.add "" @@ -125,3 +133,66 @@ proc buildReport*(records: seq[Record], source: string, outLines.add "amino acid codes. A '*' marks a stop codon." outLines.join("\n") & "\n" + +proc buildFastqReport*(records: seq[FastqRecord], source: string, + options = ReportOptions(), generatedAt = ""): string = + ## Build a text report for a list of FASTQ records. + ## + ## `source` names the origin of the data, such as a file path. When + ## `generatedAt` is empty, the current UTC time is used. + let whenGenerated = + if generatedAt.len > 0: generatedAt + else: now().utc.format("yyyy-MM-dd'T'HH:mm:ss'Z'") + + var + outLines: seq[string] + allCounts = NucleotideCounts() + + for record in records: + let counts = nucleotideCounts(record.sequence) + allCounts.a += counts.a + allCounts.c += counts.c + allCounts.g += counts.g + allCounts.t += counts.t + allCounts.u += counts.u + allCounts.n += counts.n + + outLines.add "Codon Atlas report" + outLines.add "==================" + outLines.add "Source: " & source + outLines.add "Generated: " & whenGenerated + outLines.add "Records: " & $records.len + outLines.add "Total canonical bases: " & $totalBases(allCounts) + outLines.add "" + + outLines.add "Nucleotide composition" + outLines.add "----------------------" + for index, record in records: + let counts = nucleotideCounts(record.sequence) + outLines.add "Record " & $(index + 1) & ": " & record.id + if record.description.len > 0: + outLines.add " Description: " & record.description + outLines.add " Length: " & $totalBases(counts) & " canonical bases" + outLines.add " Bases: " & formatCounts(counts) + outLines.add " GC content: " & fmtDecimal(gcPercent(record.sequence), 1) & "%" + let summary = qualitySummary(record.quality) + outLines.add " Quality: mean " & fmtDecimal(summary.mean, 1) & + ", min " & $summary.min & ", max " & $summary.max + if options.qualityThreshold > 0: + let below = basesBelowQuality(record.quality, options.qualityThreshold) + outLines.add " Bases below Q" & $options.qualityThreshold & ": " & $below + outLines.add "" + + outLines.add "Overall" + outLines.add "-------" + outLines.add "Total canonical bases: " & $totalBases(allCounts) + outLines.add "GC content: " & fmtDecimal(gcPercentOf(allCounts), 1) & "%" + + if options.kmer > 0: + outLines.add kmerSection(records.mapIt(it.sequence), options.kmer, options.kmerTop) + + outLines.add "" + outLines.add "Quality scores use the Phred+33 encoding. A higher score" + outLines.add "means a more confident base call." + + outLines.join("\n") & "\n" diff --git a/src/codonatlas_cli.nim b/src/codonatlas_cli.nim index 4b35b82..4447c4f 100644 --- a/src/codonatlas_cli.nim +++ b/src/codonatlas_cli.nim @@ -7,32 +7,42 @@ import std/strutils import codonatlas const - Version = "0.1.0" + Version = "0.2.0" ToolName = "codonatlas_cli" proc printUsage() = echo ToolName & " " & Version & " - composition statistics for sequence files" echo "" echo "Usage:" - echo " " & ToolName & " [options]" + echo " " & ToolName & " [options]" echo "" echo "Options:" - echo " -o, --out= Write the report to a file." - echo " -k, --kmer= Add a k-mer count section for k = n." - echo " -t, --top= Show only the top n k-mers." - echo " -h, --help Show this help." - echo " -v, --version Show the version." + echo " -o, --out= Write the report to a file." + echo " -k, --kmer= Add a k-mer count section for k = n." + echo " -t, --top= Show only the top n k-mers." + echo " -q, --quality= Show bases below Q in a FASTQ report." + echo " -f, --format= Input format: auto, fasta, or fastq." + echo " -h, --help Show this help." + echo " -v, --version Show the version." proc fail(message: string) = stderr.writeLine("Error: " & message) quit(1) +proc detectFormat(path: string): string = + let ext = toLowerAscii(path.splitFile.ext) + case ext + of ".fastq", ".fq": "fastq" + else: "fasta" + proc main() = var input: seq[string] outFile = "" + format = "auto" kmer = 0 top = 0 + quality = 0 showHelp = false showVersion = false @@ -55,6 +65,13 @@ proc main() = top = parseInt(val) except ValueError: fail("invalid top value: " & val) + of "quality", "q": + try: + quality = parseInt(val) + except ValueError: + fail("invalid quality threshold: " & val) + of "format", "f": + format = toLowerAscii(val) of "help", "h": showHelp = true of "version", "v": @@ -71,22 +88,34 @@ proc main() = echo Version quit(0) if input.len != 1: - stderr.writeLine("Error: provide exactly one FASTA file.") + stderr.writeLine("Error: provide exactly one sequence file.") printUsage() quit(1) + if format notin ["auto", "fasta", "fastq"]: + fail("invalid format: " & format & " (choose auto, fasta, or fastq)") let inFile = input[0] if not fileExists(inFile): fail("file not found: " & inFile) - var records: seq[Record] - try: - records = readFasta(inFile) - except FastaError as exc: - fail("invalid FASTA in " & inFile & ": " & exc.msg) + let fmt = if format == "auto": detectFormat(inFile) else: format + let options = ReportOptions(kmer: kmer, kmerTop: top, qualityThreshold: quality) + var report: string - let options = ReportOptions(kmer: kmer, kmerTop: top) - let report = buildReport(records, source = inFile, options = options) + if fmt == "fastq": + var records: seq[FastqRecord] + try: + records = readFastq(inFile) + except FastqError as exc: + fail("invalid FASTQ in " & inFile & ": " & exc.msg) + report = buildFastqReport(records, source = inFile, options = options) + else: + var records: seq[Record] + try: + records = readFasta(inFile) + except FastaError as exc: + fail("invalid FASTA in " & inFile & ": " & exc.msg) + report = buildReport(records, source = inFile, options = options) if outFile.len > 0: writeFile(outFile, report) diff --git a/tests/data/expected_fastq_report.txt b/tests/data/expected_fastq_report.txt new file mode 100644 index 0000000..9bd0d71 --- /dev/null +++ b/tests/data/expected_fastq_report.txt @@ -0,0 +1,39 @@ +Codon Atlas report +================== +Source: sample.fastq +Generated: 2026-08-02T12:00:00Z +Records: 2 +Total canonical bases: 35 + +Nucleotide composition +---------------------- +Record 1: read1 + Description: synthetic high-quality read, 23 bp + Length: 23 canonical bases + Bases: A 5 C 5 G 7 T 6 U 0 N 0 + GC content: 52.2% + Quality: mean 38.3, min 37, max 40 + Bases below Q20: 0 + +Record 2: read2 + Description: synthetic low-quality read, 12 bp + Length: 12 canonical bases + Bases: A 3 C 4 G 2 T 3 U 0 N 0 + GC content: 50.0% + Quality: mean 7.0, min 0, max 13 + Bases below Q20: 12 + +Overall +------- +Total canonical bases: 35 +GC content: 51.4% + +K-mer counts (k = 3) +-------------------------------- +Top 3 of 25 distinct k-mers. +ATG 2 +CCG 2 +CGT 2 + +Quality scores use the Phred+33 encoding. A higher score +means a more confident base call. diff --git a/tests/data/sample.fastq b/tests/data/sample.fastq new file mode 100644 index 0000000..9ef65e4 --- /dev/null +++ b/tests/data/sample.fastq @@ -0,0 +1,8 @@ +@read1 synthetic high-quality read, 23 bp +ATGCGTACGTTAGCCGATGCAGT ++ +IIIIIHHHHHHHHFFFFFFFFFF +@read2 synthetic low-quality read, 12 bp +AATTCCGGATCC ++ +!#$%&()*+,-. diff --git a/tests/test_fastq.nim b/tests/test_fastq.nim new file mode 100644 index 0000000..7376e04 --- /dev/null +++ b/tests/test_fastq.nim @@ -0,0 +1,166 @@ +import std/os +import std/strutils +import std/unittest + +import codonatlas + +proc fixturePath(): string = + currentSourcePath().parentDir() / "data" / "sample.fastq" + +proc goldenPath(): string = + currentSourcePath().parentDir() / "data" / "expected_fastq_report.txt" + +suite "quality scores": + test "Phred+33 encoding": + check qualityScore('!') == 0 + check qualityScore('I') == 40 + check qualityScore('~') == 93 + check qualityScore('@', offset = 64) == 0 + + test "out-of-range scores are rejected": + expect ValueError: + discard qualityScore(' ') + +suite "quality summary": + test "known strings": + let summary = qualitySummary("IHH") + check summary.count == 3 + check summary.min == 39 + check summary.max == 40 + check abs(summary.mean - 39.3333) < 0.001 + + test "empty string yields zeros": + let summary = qualitySummary("") + check summary.count == 0 + check summary.min == 0 + check summary.max == 0 + check summary.mean == 0.0 + +suite "bases below threshold": + test "strictly below the threshold": + check basesBelowQuality("IHH", 40) == 2 + check basesBelowQuality("IHH", 39) == 0 + check basesBelowQuality("", 20) == 0 + +suite "FASTQ parsing": + test "multiple records with descriptions": + let records = parseFastq("@a\nAT\n+\nII\n@b note here\nGC\n+\n!!\n") + check records.len == 2 + check records[0].id == "a" + check records[0].description == "" + check records[0].sequence == "AT" + check records[0].quality == "II" + check records[1].id == "b" + check records[1].description == "note here" + check records[1].sequence == "GC" + check records[1].quality == "!!" + + test "wrapped sequence and quality lines": + let records = parseFastq("@a\nAT\nGC\n+\nII\n!!\n") + check records[0].sequence == "ATGC" + check records[0].quality == "II!!" + + test "blank lines between records are ignored": + let records = parseFastq("@a\nAT\n+\nII\n\n@b\nGC\n+\n!!\n") + check records.len == 2 + + test "whitespace inside lines is removed": + let records = parseFastq("@a\nA T\nGC\n+\nI I\n!!\n") + check records[0].sequence == "ATGC" + check records[0].quality == "II!!" + + test "lowercase is preserved": + let records = parseFastq("@a\natgc\n+\nIIII\n") + check records[0].sequence == "atgc" + + test "the plus line may carry text": + let records = parseFastq("@a\nAT\n+ some text\nII\n") + check records[0].sequence == "AT" + check records[0].quality == "II" + + test "empty input produces no records": + check parseFastq("").len == 0 + + test "sequence before a header is rejected": + expect FastqError: + discard parseFastq("ATGC\n+\nIIII\n") + + test "empty header is rejected": + expect FastqError: + discard parseFastq("@\nAT\n+\nII\n") + + test "record without sequence is rejected": + expect FastqError: + discard parseFastq("@a\n+\nII\n") + + test "missing plus separator is rejected": + expect FastqError: + discard parseFastq("@a\nAT\nII\n") + + test "missing quality line is rejected": + expect FastqError: + discard parseFastq("@a\nAT\n+\n") + + test "quality shorter than sequence is rejected": + expect FastqError: + discard parseFastq("@a\nAT\n+\nI\n") + + test "quality longer than sequence is rejected": + expect FastqError: + discard parseFastq("@a\nAT\n+\nIII\n") + + test "readFastq reads the bundled fixture": + let records = readFastq(fixturePath()) + check records.len == 2 + check records[0].id == "read1" + check records[0].sequence.len == 23 + check records[0].quality.len == 23 + check records[1].id == "read2" + check records[1].description.len > 0 + +suite "FASTQ report": + test "required sections are present": + let records = readFastq(fixturePath()) + let report = buildFastqReport(records, source = "sample.fastq", + generatedAt = "2026-08-02T12:00:00Z") + check report.contains("Codon Atlas report") + check report.contains("Source: sample.fastq") + check report.contains("Generated: 2026-08-02T12:00:00Z") + check report.contains("Records: 2") + check report.contains("Total canonical bases: 35") + check report.contains("Quality: mean 38.3, min 37, max 40") + + test "the below-threshold line appears on request": + let records = readFastq(fixturePath()) + let report = buildFastqReport(records, source = "sample.fastq", + options = ReportOptions(qualityThreshold: 20), + generatedAt = "2026-08-02T12:00:00Z") + check report.contains("Bases below Q20: 12") + + test "the below-threshold line stays hidden by default": + let records = readFastq(fixturePath()) + let report = buildFastqReport(records, source = "sample.fastq", + generatedAt = "2026-08-02T12:00:00Z") + check not report.contains("Bases below Q20") + + test "k-mer section appears on request": + let records = readFastq(fixturePath()) + let report = buildFastqReport(records, source = "sample.fastq", + options = ReportOptions(kmer: 3), + generatedAt = "2026-08-02T12:00:00Z") + check report.contains("K-mer counts (k = 3)") + + test "output is deterministic": + let records = readFastq(fixturePath()) + let first = buildFastqReport(records, source = "sample.fastq", + generatedAt = "2026-08-02T12:00:00Z") + let second = buildFastqReport(records, source = "sample.fastq", + generatedAt = "2026-08-02T12:00:00Z") + check first == second + + test "golden report matches the bundled fixture": + let records = readFastq(fixturePath()) + let report = buildFastqReport(records, source = "sample.fastq", + options = ReportOptions(kmer: 3, kmerTop: 3, qualityThreshold: 20), + generatedAt = "2026-08-02T12:00:00Z") + check report == readFile(goldenPath())