From ddf8178a8dd6778087bc2fd9c7ec53109f05849d Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 21 Oct 2026 18:34:51 +0000 Subject: [PATCH 1/6] fix crash/dead-flag bugs in 10X cmd wiring The "add 10X support for basic FQ functions" commit registered several new flags on the wrong command, which panics the entire binary on startup regardless of which subcommand is invoked: - cmd/count.go and cmd/extract.go each registered their new --bc-list flag on convertFqCmd instead of their own Fq command, which also means convertFqCmd tried to register "bc-list" three times total (once correctly in convert.go, twice more here) -- pflag panics on a redefined flag. - cmd/stats.go registered --bc-list/--threads on unlinkedFqCmd instead of statsFqCmd, on top of rm_unlinked.go's own (correct) registration of those same flags on unlinkedFqCmd -- another redefinition panic, and statsFqCmd ends up with no --bc-list/--threads flags at all. - cmd/rm_unlinked.go gave unlinkedFqCmd two flags under the same -b shorthand (--bc-count, a bool from an earlier change, and the new --bc-list) -- pflag panics on a duplicate shorthand. Its Args func also still called GetString("bc-count") after that flag became a bool, which would have been the next failure once the panic was fixed. Also added the missing "10x"/"10X" bc-list sentinel exemption to the Args validators that skip filecheck for a literal whitelist path (cmd/sort.go, cmd/ncbi.go, cmd/count.go, cmd/stats.go) -- without it, `--bc-list 10x` (meaning "use the embedded official whitelist") was rejected before ever reaching fastq.Load10X, which already understood that sentinel. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_015tF7ppXiNo2Adbv8YXcwGL --- cmd/count.go | 13 ++++++++++++- cmd/extract.go | 2 +- cmd/rm_unlinked.go | 13 +------------ cmd/stats.go | 22 +++++++++++++++++----- 4 files changed, 31 insertions(+), 19 deletions(-) diff --git a/cmd/count.go b/cmd/count.go index 0d0d6e5..6a7901d 100644 --- a/cmd/count.go +++ b/cmd/count.go @@ -63,6 +63,17 @@ var countFqCmd = &cobra.Command{ return err } } + mapfile, err := cmd.Flags().GetString("bc-list") + if err != nil { + return err + } + switch mapfile { + case "", "10x", "10X": // skip check + default: + if err := filecheck(mapfile); err != nil { + return err + } + } return nil }, RunE: func(cmd *cobra.Command, args []string) error { @@ -83,7 +94,7 @@ func init() { countXamCmd.Flags().BoolP("invalid", "i", false, "Include invalid barcodes") countXamCmd.Flags().IntP("threads", "@", 2, "Threads to use") - convertFqCmd.Flags().StringP("bc-list", "b", "", "Valid barcode list [if 10X input]") + countFqCmd.Flags().StringP("bc-list", "b", "", "Valid barcode list [if 10X input]") countFqCmd.Flags().BoolP("invalid", "i", false, "Include invalid barcodes") countFqCmd.Flags().IntP("threads", "@", 2, "Threads to use") } diff --git a/cmd/extract.go b/cmd/extract.go index 9b25bff..cddf353 100644 --- a/cmd/extract.go +++ b/cmd/extract.go @@ -96,7 +96,7 @@ func init() { extractXamCmd.Flags().BoolP("invalid", "i", false, "Include invalid barcodes") extractXamCmd.Flags().IntP("threads", "@", 2, "Threads to use") - convertFqCmd.Flags().StringP("bc-list", "b", "", "Valid barcode list [if 10X input]") + extractFqCmd.Flags().StringP("bc-list", "b", "", "Valid barcode list [if 10X input]") extractFqCmd.Flags().BoolP("invalid", "i", false, "Include invalid barcodes") extractFqCmd.Flags().IntP("threads", "@", 2, "Threads to use") } diff --git a/cmd/rm_unlinked.go b/cmd/rm_unlinked.go index 38a82f3..faccc71 100644 --- a/cmd/rm_unlinked.go +++ b/cmd/rm_unlinked.go @@ -101,17 +101,6 @@ var unlinkedFqCmd = &cobra.Command{ return err } } - bcCount, err := cmd.Flags().GetString("bc-count") - if err != nil { - return err - } - if bcCount != "" { - err = ensureWritableDir(bcCount) - if err != nil { - return err - } - } - mapfile, err := cmd.Flags().GetString("bc-list") if err != nil { return err @@ -149,7 +138,7 @@ func init() { //TODO SINGLETONS SHOULD BE A FLAG, OUTPUT SHOULD BE PREFIX.SINGLE.... unlinkedFqCmd.Flags().BoolP("unlinked", "u", false, "Write valid singleton records to prefix.unlinked (optional)") - unlinkedFqCmd.Flags().BoolP("bc-count", "b", false, "Write valid barcode counts to prefix.counts (optional)") + unlinkedFqCmd.Flags().Bool("bc-count", false, "Write valid barcode counts to prefix.counts (optional)") unlinkedFqCmd.Flags().StringP("bc-list", "b", "", "Valid barcode list [if 10X input]") unlinkedFqCmd.Flags().IntP("threads", "@", 2, "Threads to use") } diff --git a/cmd/stats.go b/cmd/stats.go index d55869d..9a61e62 100644 --- a/cmd/stats.go +++ b/cmd/stats.go @@ -89,9 +89,21 @@ var statsFqCmd = &cobra.Command{ } switch strings.ToLower(platform) { case "", stats.Haplotagging, stats.Stlfr, stats.Tellseq: - return nil + default: + return fmt.Errorf("unknown platform %q, expected haplotagging, stlfr or tellseq", platform) } - return fmt.Errorf("unknown platform %q, expected haplotagging, stlfr or tellseq", platform) + mapfile, err := cmd.Flags().GetString("bc-list") + if err != nil { + return err + } + switch mapfile { + case "", "10x", "10X": // skip check + default: + if err := filecheck(mapfile); err != nil { + return err + } + } + return nil }, RunE: func(cmd *cobra.Command, args []string) error { platform, _ := cmd.Flags().GetString("platform") @@ -111,7 +123,7 @@ func init() { statsXamCmd.Flags().IntP("distance-threshold", "d", 50000, "Distance (bp) that splits a barcode's reads into separate molecules, where 0 skips deconvolution") statsXamCmd.Flags().IntP("threads", "@", 2, "Threads to use") - statsFqCmd.Flags().StringP("platform", "p", "", "Linked-read technology: haplotagging, stlfr or tellseq (default: detect from the first 100 reads)") - unlinkedFqCmd.Flags().StringP("bc-list", "b", "", "Valid barcode list [if 10X input]") - unlinkedFqCmd.Flags().IntP("threads", "@", 2, "Threads to use") + statsFqCmd.Flags().StringP("platform", "p", "", "Linked-read technology: haplotagging, stlfr or tellseq (default: detect from the first 100 reads). For 10X data, use --bc-list instead.") + statsFqCmd.Flags().StringP("bc-list", "b", "", "Valid barcode list [if 10X input]") + statsFqCmd.Flags().IntP("threads", "@", 2, "Threads to use") } From 49d14e15596f0ec092841308d7ba86ad60e7582b Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 21 Oct 2026 18:35:03 +0000 Subject: [PATCH 2/6] fix remaining 10X sentinel/doc gaps in cmd and fastq - cmd/ncbi.go, cmd/sort.go: same missing "10x"/"10X" sentinel exemption as the previous commit, for the two FASTQ commands whose --bc-list support predates this feature's addition elsewhere. - cmd/root.go: fqCmd's Long text still named the flag --bc-map, which was renamed to --bc-list; fixed the name and a typo ("offical"). - fastq/parser_10x.go: typo in Load10X's doc comment ("rease" the bcmap file). Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_015tF7ppXiNo2Adbv8YXcwGL --- cmd/ncbi.go | 8 ++++++-- cmd/root.go | 2 +- cmd/sort.go | 4 +++- fastq/parser_10x.go | 2 +- 4 files changed, 11 insertions(+), 5 deletions(-) diff --git a/cmd/ncbi.go b/cmd/ncbi.go index 82ec4a6..f6c24ca 100644 --- a/cmd/ncbi.go +++ b/cmd/ncbi.go @@ -69,8 +69,12 @@ var ncbiFq2XamCmd = &cobra.Command{ if err != nil { return err } - if mapfile != "" && filecheck(mapfile) != nil { - return err + switch mapfile { + case "", "10x", "10X": // skip check + default: + if filecheck(mapfile) != nil { + return err + } } return nil }, diff --git a/cmd/root.go b/cmd/root.go index 7407045..ac2bfdd 100644 --- a/cmd/root.go +++ b/cmd/root.go @@ -26,7 +26,7 @@ var samCmd = &cobra.Command{ var fqCmd = &cobra.Command{ Use: "fq", Short: "FASTQ operations", - Long: "Input can be one FASTQ file or file pair. If input is 10X data, provide a barcode 'whitelist' to --bc-map or use the embedded offical 10X whitelist with '--bc-map 10x'.", + Long: "Input can be one FASTQ file or file pair. If input is 10X data, provide a barcode whitelist to --bc-list, or use the embedded official 10X whitelist with '--bc-list 10x'.", } func init() { diff --git a/cmd/sort.go b/cmd/sort.go index 23c90c4..d211030 100644 --- a/cmd/sort.go +++ b/cmd/sort.go @@ -81,7 +81,9 @@ var sortFqCmd = &cobra.Command{ if err != nil { return err } - if bcMap != "" { + switch bcMap { + case "", "10x", "10X": // skip check + default: if err := filecheck(bcMap); err != nil { return err } diff --git a/fastq/parser_10x.go b/fastq/parser_10x.go index f3201f7..617678a 100644 --- a/fastq/parser_10x.go +++ b/fastq/parser_10x.go @@ -60,7 +60,7 @@ func TenXR22Corefq(rec *fastx.Record) (core CoreFq) { } // Load a 10X barcode whitelist. If bcmap is '10x' or '10X', loads the official list -// embedded in djinn, otherwise rease the bcmap file +// embedded in djinn, otherwise reads the bcmap file func Load10X(bcmap string) (barcodeMap map[string]struct{}, err error) { switch bcmap { case "10X", "10x": From 417fccd11aad4e0965969c83edfddefe04e602cb Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 21 Oct 2026 18:35:17 +0000 Subject: [PATCH 3/6] fix three broken/missing 10X code paths - stats/stats.fq.go: StatsFQ built a 10X parser/scheme from --bc-list into parse/sch, then unconditionally called "parse, sch, err := chooseParser(fq, platform)" right after -- chooseParser only knows haplotagging/stlfr/tellseq, and reusing parse/sch here (same block as their earlier `var` decls, so this isn't shadowing) discarded the 10X values every time regardless of --bc-list. The entire 10X branch was dead code; `djinn fq stats` always fell through to auto-detection, which 10X data can't satisfy (its barcode is just an unmarked sequence prefix). - invalid/invalid.fq.go: FilterInvalid10XFQ's paired branch never wrote a valid R2 record to R2Writer at all -- the only call to r2core.ToTenX(R2Writer.Writer) was inside the "!r1core.Valid" branch, so prefix.R2.fq.gz came out empty for every properly-paired, valid-barcode 10X input. It also had no separate invalid-output file for 10X (unlike every other format): "kept" invalid records were written into the same file as valid ones instead of their own prefix.invalid.R*.fq.gz. Added R1WriterInv/R2WriterInv and fixed the R2 write to happen unconditionally for the valid case. - unlinked/unlinked.fq.go: RemoveUnlinkedFQ accepted a bcmap parameter (added so rm-unlinked's new --bc-list flag would have somewhere to go) but never used it anywhere in the function body -- 10X support was wired up at the CLI layer and silently absent underneath it. Added RemoveUnlinked10XFQ (and its own getFqCount10X, since a 10X barcode lives only in R1) mirroring the single-end/paired structure already used by FilterInvalid10XFQ/ConvertFrom10X, and dispatch to it from RemoveUnlinkedFQ when bcmap != "". - extract/extract.fq.go: Extract10XFQ duplicated the whitelist-sentinel switch that fastq.Load10X already implements (and that this same "add 10X support" commit introduced), instead of calling it like Count10XFQ/ConvertFrom10X/FilterInvalid10XFQ/StatsFQ all do. No functional difference, just inconsistent with its own siblings. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_015tF7ppXiNo2Adbv8YXcwGL --- extract/extract.fq.go | 19 +---- invalid/invalid.fq.go | 42 ++++++++-- stats/stats.fq.go | 12 +-- unlinked/unlinked.fq.go | 178 ++++++++++++++++++++++++++++++++++++++++ 4 files changed, 221 insertions(+), 30 deletions(-) diff --git a/extract/extract.fq.go b/extract/extract.fq.go index c25639d..35a9647 100644 --- a/extract/extract.fq.go +++ b/extract/extract.fq.go @@ -6,7 +6,6 @@ import ( "io" "os" - "github.com/blinkseq/djinn/barcodes" "github.com/blinkseq/djinn/fastq" "github.com/shenwei356/bio/seq" @@ -77,23 +76,11 @@ func ExtractFQ(fqs []string, invalid bool, bcmap string) error { func Extract10XFQ(fqs []string, invalid bool, bcmap string) error { seq.ValidateSeq = false - var bclist map[string]struct{} - var err error - - switch bcmap { - case "10X", "10x": - bclist, err = barcodes.LoadTenXSet() - if err != nil { - return err - } - default: - bclist, err = barcodes.LoadTenXFromFile(bcmap) - if err != nil { - return err - } + bclist, err := fastq.Load10X(bcmap) + if err != nil { + return err } - seq.ValidateSeq = false // ── open barcode writer ─────────────────────────────────────────────────────────── writer := bufio.NewWriter(os.Stdout) defer writer.Flush() diff --git a/invalid/invalid.fq.go b/invalid/invalid.fq.go index ae43024..db56ab5 100644 --- a/invalid/invalid.fq.go +++ b/invalid/invalid.fq.go @@ -112,6 +112,17 @@ func FilterInvalid10XFQ(fqs []string, prefix string, keepInvalid bool, bcmap str } defer R1Writer.Close() + // Invalid records go to their own prefix.invalid.R*.fq.gz files, same + // as the non-10X path -- not mixed back into the main output. + var R1WriterInv *xopen.Writer + if keepInvalid { + R1WriterInv, err = xopen.Wopen(prefix + ".invalid.R1.fq.gz") + if err != nil { + return err + } + defer R1WriterInv.Close() + } + if len(fqs) == 1 { // Single-end 10X input: no R2 to propagate the barcode onto. for { @@ -125,7 +136,7 @@ func FilterInvalid10XFQ(fqs []string, prefix string, keepInvalid bool, bcmap str r1core, ok := fastq.TenXR12Corefq(_r1, bclist) if !ok || !r1core.Valid { if keepInvalid { - r1core.ToTenX(R1Writer.Writer) + r1core.ToTenX(R1WriterInv.Writer) } continue } @@ -144,6 +155,15 @@ func FilterInvalid10XFQ(fqs []string, prefix string, keepInvalid bool, bcmap str } defer R2Writer.Close() + var R2WriterInv *xopen.Writer + if keepInvalid { + R2WriterInv, err = xopen.Wopen(prefix + ".invalid.R2.fq.gz") + if err != nil { + return err + } + defer R2WriterInv.Close() + } + var r1core, r2core fastq.CoreFq Loop: @@ -158,7 +178,7 @@ func FilterInvalid10XFQ(fqs []string, prefix string, keepInvalid bool, bcmap str r1core, ok = fastq.TenXR12Corefq(_r1, bclist) if !ok || !r1core.Valid { if keepInvalid { - r1core.ToTenX(R1Writer.Writer) + r1core.ToTenX(R1WriterInv.Writer) } continue } @@ -177,19 +197,25 @@ func FilterInvalid10XFQ(fqs []string, prefix string, keepInvalid bool, bcmap str } case nil: r2core = fastq.TenXR22Corefq(_r2) - // check if R1 exists and match the seq ID + // check if R1 exists and match the seq ID; a record with + // no matching R1 (should only happen if the files are + // mismatched/truncated) has no barcode info and is + // therefore treated the same as an invalid barcode. + valid := false if r1err == nil { if ok := slices.Equal(r1core.ID, r2core.ID); !ok { return fmt.Errorf("Misaligned R1 and R2 records. The 10X format only stores the barcode at the beginning of read 1, meaning Read 2 must have the same sequence ID as Read 1 to reliably assign a barcode to Read 2, but they do not match.\nRead 1:\n%s\nRead 2:\n%s", r1core.ID, r2core.ID) } r2core.Barcode = r1core.Barcode - if !r1core.Valid { - if keepInvalid { - r2core.ToTenX(R2Writer.Writer) - } - continue + valid = r1core.Valid + } + if !valid { + if keepInvalid { + r2core.ToTenX(R2WriterInv.Writer) } + continue } + r2core.ToTenX(R2Writer.Writer) default: return r2err } diff --git a/stats/stats.fq.go b/stats/stats.fq.go index d5e6092..5bbe94b 100644 --- a/stats/stats.fq.go +++ b/stats/stats.fq.go @@ -43,18 +43,18 @@ func StatsFQ(fq, platform string, bcMap string) error { var sch scheme if bcMap != "" { - var bclist map[string]struct{} bclist, err := fastq.Load10X(bcMap) if err != nil { return err } parse = func(rec *fastx.Record) (fastq.CoreFq, bool) { return fastq.TenXR12Corefq(rec, bclist) } sch = tenxScheme - } - - parse, sch, err := chooseParser(fq, platform) - if err != nil { - return err + } else { + var err error + parse, sch, err = chooseParser(fq, platform) + if err != nil { + return err + } } seq.ValidateSeq = false // reads can hold N and other ambiguity codes; nothing here needs them checked diff --git a/unlinked/unlinked.fq.go b/unlinked/unlinked.fq.go index 534ff90..ac6bbf9 100644 --- a/unlinked/unlinked.fq.go +++ b/unlinked/unlinked.fq.go @@ -6,6 +6,7 @@ import ( "io" "maps" "os" + "slices" "strconv" "github.com/blinkseq/djinn/fastq" @@ -86,7 +87,184 @@ func getFqCount(infiles []string) (map[string]int, error) { return set, nil } +// getFqCount10X counts valid 10X barcodes from R1 only -- unlike the other +// linked-read formats, a 10X barcode lives solely in R1's sequence prefix, +// so R2 has nothing to independently count. +func getFqCount10X(fqs []string, bclist map[string]struct{}) (map[string]int, error) { + seq.ValidateSeq = false + set := make(map[string]int, 4_000_000) + + fqReader, err := fastx.NewReader(seq.DNA, fqs[0], "") + if err != nil { + return nil, fmt.Errorf("opening %s: %w", fqs[0], err) + } + defer fqReader.Close() + + for { + rec, err := fqReader.Read() + if err == io.EOF { + break + } + if err != nil { + return nil, err + } + core, ok := fastq.TenXR12Corefq(rec, bclist) + if !ok || !core.Valid { + continue + } + set[string(core.Barcode)]++ + } + return set, nil +} + +// RemoveUnlinked10XFQ is the 10X-input counterpart of RemoveUnlinkedFQ: a +// barcode is "linked" if its R1 sequence-prefix matches the whitelist on 2 +// or more reads, same threshold as every other format. +func RemoveUnlinked10XFQ(fqs []string, prefix string, keepSingle, barcodecount bool, bcmap string) error { + bclist, err := fastq.Load10X(bcmap) + if err != nil { + return err + } + + bcCounts, err := getFqCount10X(fqs, bclist) + if err != nil { + return err + } + if barcodecount { + f, err := os.Create(prefix + ".count") + if err != nil { + return err + } + writer := bufio.NewWriter(f) + for key, val := range bcCounts { + writer.WriteString(key) + writer.WriteByte('\t') + writer.WriteString(strconv.Itoa(val)) + writer.WriteByte('\n') + } + writer.Flush() + f.Close() + } + + R1Reader, err := fastx.NewReader(seq.DNA, fqs[0], "") + if err != nil { + return fmt.Errorf("opening %s: %w", fqs[0], err) + } + defer R1Reader.Close() + + R1Writer, err := xopen.Wopen(prefix + ".R1.fq.gz") + if err != nil { + return err + } + defer R1Writer.Close() + + var R1WriterSingle *xopen.Writer + if keepSingle { + R1WriterSingle, err = xopen.Wopen(prefix + ".unlinked.R1.fq.gz") + if err != nil { + return err + } + defer R1WriterSingle.Close() + } + + route := func(core *fastq.CoreFq, w, wSingle *xopen.Writer) { + switch { + case bcCounts[string(core.Barcode)] >= 2: + core.ToTenX(w.Writer) + case keepSingle: + core.ToTenX(wSingle.Writer) + } + } + + if len(fqs) == 1 { + // Single-end 10X input: no R2 to propagate the barcode onto. + for { + _r1, err := R1Reader.Read() + if err == io.EOF { + break + } + if err != nil { + return err + } + r1core, ok := fastq.TenXR12Corefq(_r1, bclist) + if !ok || !r1core.Valid { + continue + } + route(&r1core, R1Writer, R1WriterSingle) + } + return nil + } + + R2Reader, err := fastx.NewReader(seq.DNA, fqs[1], "") + if err != nil { + return fmt.Errorf("opening %s: %w", fqs[1], err) + } + defer R2Reader.Close() + + R2Writer, err := xopen.Wopen(prefix + ".R2.fq.gz") + if err != nil { + return err + } + defer R2Writer.Close() + + var R2WriterSingle *xopen.Writer + if keepSingle { + R2WriterSingle, err = xopen.Wopen(prefix + ".unlinked.R2.fq.gz") + if err != nil { + return err + } + defer R2WriterSingle.Close() + } + + var r1core, r2core fastq.CoreFq + +Loop: + for { + // process R1, which doesn't need an R2 + _r1, r1err := R1Reader.Read() + switch r1err { + case io.EOF: + + case nil: + var ok bool + r1core, ok = fastq.TenXR12Corefq(_r1, bclist) + if !ok || !r1core.Valid { + r1core.Barcode = nil // marks this pair as having no usable barcode for R2 below + } else { + route(&r1core, R1Writer, R1WriterSingle) + } + default: + return r1err + } + + // process R2 + _r2, r2err := R2Reader.Read() + switch r2err { + case io.EOF: + // break when both are exhausted + if r1err == io.EOF { + break Loop + } + case nil: + r2core = fastq.TenXR22Corefq(_r2) + if r1err == nil && len(r1core.Barcode) > 0 { + if ok := slices.Equal(r1core.ID, r2core.ID); !ok { + return fmt.Errorf("Misaligned R1 and R2 records. The 10X format only stores the barcode at the beginning of read 1, meaning Read 2 must have the same sequence ID as Read 1 to reliably assign a barcode to Read 2, but they do not match.\nRead 1:\n%s\nRead 2:\n%s", r1core.ID, r2core.ID) + } + r2core.Barcode = r1core.Barcode + route(&r2core, R2Writer, R2WriterSingle) + } + default: + return r2err + } + } + return nil +} + func RemoveUnlinkedFQ(fqs []string, prefix string, keepSingle, barcodecount bool, bcmap string) error { + if bcmap != "" { + return RemoveUnlinked10XFQ(fqs, prefix, keepSingle, barcodecount, bcmap) + } // guard against draining stdin when getting barcode counts bcCounts, err := getFqCount(fqs) if err != nil { From a7ed8aab505a5518baaf17808f0f45687e474dad Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 21 Oct 2026 18:35:26 +0000 Subject: [PATCH 4/6] tests: fix signature drift, add 10X regression coverage CountFQ, ExtractFQ, FilterInvalidFQ, StatsFQ, and RemoveUnlinkedFQ all gained a new parameter in the "add 10X support" commit; every existing call site in these five packages' test files was still using the old signature, so none of them even compiled. None of count/, extract/, invalid/, or unlinked/ had any 10X test coverage at all before this, which is exactly how the bugs fixed in the previous commit went unnoticed. Added: - count: Count10XFQ end-to-end (dedup + unrecognized-barcode drop), plus the "10x" embedded-whitelist sentinel. - extract: Extract10XFQ end-to-end, plus confirms a second (R2) file is correctly ignored rather than misread. - invalid: regression tests for both FilterInvalid10XFQ bugs just fixed (an empty R2 output, and invalid records leaking into the main output instead of their own file), plus single-end. - unlinked: RemoveUnlinked10XFQ singleton/linked routing (paired and single-end) -- this function didn't exist before this branch. - stats: a regression test for the dead-10X-branch bug just fixed, asserting real Barcodes_Valid output rather than just "no error". Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_015tF7ppXiNo2Adbv8YXcwGL --- count/count_fq_test.go | 87 ++++++++++++++++++++-- extract/extract_fq_test.go | 79 +++++++++++++++++++- invalid/invalid_fq_test.go | 139 ++++++++++++++++++++++++++++++++++- stats/stats_fq_test.go | 42 ++++++++++- unlinked/unlinked_fq_test.go | 111 ++++++++++++++++++++++++++-- 5 files changed, 436 insertions(+), 22 deletions(-) diff --git a/count/count_fq_test.go b/count/count_fq_test.go index e761175..4fbf110 100644 --- a/count/count_fq_test.go +++ b/count/count_fq_test.go @@ -1,7 +1,9 @@ package count import ( + "os" "path/filepath" + "strings" "testing" ) @@ -22,7 +24,7 @@ func TestCountFQ_R1Only(t *testing.T) { writeHapFastq(t, r1, reads) out, err := captureStdout(t, func() error { - return CountFQ([]string{r1}, false) + return CountFQ([]string{r1}, false, "") }) if err != nil { t.Fatalf("CountFQ returned error: %v", err) @@ -52,7 +54,7 @@ func TestCountFQ_IncludeInvalid(t *testing.T) { writeHapFastq(t, r1, reads) out, err := captureStdout(t, func() error { - return CountFQ([]string{r1}, true) + return CountFQ([]string{r1}, true, "") }) if err != nil { t.Fatalf("CountFQ returned error: %v", err) @@ -84,7 +86,7 @@ func TestCountFQ_R1R2Union(t *testing.T) { writeHapFastq(t, r2, r2Reads) out, err := captureStdout(t, func() error { - return CountFQ([]string{r1, r2}, false) + return CountFQ([]string{r1, r2}, false, "") }) if err != nil { t.Fatalf("CountFQ returned error: %v", err) @@ -109,7 +111,7 @@ func TestCountFQ_FewerThan100Records(t *testing.T) { writeHapFastq(t, r1, []hapRead{{bc: "AAAA", vx: 1}}) _, err := captureStdout(t, func() error { - return CountFQ([]string{r1}, false) + return CountFQ([]string{r1}, false, "") }) if err != nil { t.Fatalf("expected no error for a fastq file with fewer than 100 records, got %v", err) @@ -118,9 +120,84 @@ func TestCountFQ_FewerThan100Records(t *testing.T) { func TestCountFQ_MissingFile(t *testing.T) { _, err := captureStdout(t, func() error { - return CountFQ([]string{"/nonexistent/path/does-not-exist.fq"}, false) + return CountFQ([]string{"/nonexistent/path/does-not-exist.fq"}, false, "") }) if err == nil { t.Fatal("expected an error opening a nonexistent file, got nil") } } + +// writeTenXFastq writes a plain-text FASTQ file of 10X-format records: a +// 16bp barcode inline as the sequence prefix, no separate tag. +func writeTenXFastq(t *testing.T, path string, ids []string, barcode string) { + t.Helper() + var sb strings.Builder + for _, id := range ids { + seq := barcode + "ACGTACGTAC" + sb.WriteString("@" + id + "\n" + seq + "\n+\n" + strings.Repeat("I", len(seq)) + "\n") + } + if err := os.WriteFile(path, []byte(sb.String()), 0o644); err != nil { + t.Fatal(err) + } +} + +func writeTenXWhitelist(t *testing.T, path string, barcodes ...string) { + t.Helper() + if err := os.WriteFile(path, []byte(strings.Join(barcodes, "\n")+"\n"), 0o644); err != nil { + t.Fatal(err) + } +} + +func TestCount10XFQ(t *testing.T) { + dir := t.TempDir() + r1 := filepath.Join(dir, "r1.fq") + wl := filepath.Join(dir, "whitelist.txt") + + writeTenXWhitelist(t, wl, "AAAACCCCGGGGTTTT") + // two reads with the known barcode, one with an unrecognized prefix + writeTenXFastq(t, r1, []string{"read1", "read2"}, "AAAACCCCGGGGTTTT") + writeTenXFastq(t, filepath.Join(dir, "r1b.fq"), []string{"read3"}, "ZZZZZZZZZZZZZZZZ") + // append the unrecognized-barcode read to r1 so both are in one file + data, err := os.ReadFile(r1) + if err != nil { + t.Fatal(err) + } + extra, err := os.ReadFile(filepath.Join(dir, "r1b.fq")) + if err != nil { + t.Fatal(err) + } + if err := os.WriteFile(r1, append(data, extra...), 0o644); err != nil { + t.Fatal(err) + } + + out, err := captureStdout(t, func() error { + return CountFQ([]string{r1}, false, wl) + }) + if err != nil { + t.Fatalf("CountFQ (10X): %v", err) + } + counts := parseCounts(t, out) + if counts["AAAACCCCGGGGTTTT"] != 2 { + t.Errorf("AAAACCCCGGGGTTTT count = %d, want 2", counts["AAAACCCCGGGGTTTT"]) + } + if len(counts) != 1 { + t.Errorf("expected only the whitelisted barcode to be counted, got %v", counts) + } +} + +func TestCount10XFQ_EmbeddedWhitelistSentinel(t *testing.T) { + dir := t.TempDir() + r1 := filepath.Join(dir, "r1.fq") + writeTenXFastq(t, r1, []string{"read1"}, "AAAACCCCGGGGTTTT") + + // "10x"/"10X" means "use the embedded official whitelist" (fastq.Load10X), + // not a literal file named "10x" -- this used to be unreachable because + // the CLI's own Args validator ran filecheck("10x") before this point, + // but CountFQ/Count10XFQ themselves have always honored the sentinel. + _, err := captureStdout(t, func() error { + return CountFQ([]string{r1}, false, "10x") + }) + if err != nil { + t.Fatalf("CountFQ with embedded whitelist sentinel: %v", err) + } +} diff --git a/extract/extract_fq_test.go b/extract/extract_fq_test.go index 8fc1e0a..0d7ed36 100644 --- a/extract/extract_fq_test.go +++ b/extract/extract_fq_test.go @@ -34,7 +34,7 @@ func TestExtractFQHaplotagging(t *testing.T) { var err error out := captureStdout(t, func() { - err = ExtractFQ([]string{fq}, false) + err = ExtractFQ([]string{fq}, false, "") }) if err != nil { t.Fatalf("ExtractFQ: %v", err) @@ -59,7 +59,7 @@ func TestExtractFQIncludeInvalid(t *testing.T) { var err error out := captureStdout(t, func() { - err = ExtractFQ([]string{fq}, true) + err = ExtractFQ([]string{fq}, true, "") }) if err != nil { t.Fatalf("ExtractFQ: %v", err) @@ -84,7 +84,7 @@ func TestExtractFQMultipleFiles(t *testing.T) { var err error out := captureStdout(t, func() { - err = ExtractFQ([]string{fq1, fq2}, false) + err = ExtractFQ([]string{fq1, fq2}, false, "") }) if err != nil { t.Fatalf("ExtractFQ: %v", err) @@ -96,8 +96,79 @@ func TestExtractFQMultipleFiles(t *testing.T) { } func TestExtractFQBadFirstFile(t *testing.T) { - err := ExtractFQ([]string{"/nonexistent/path/reads.fq"}, false) + err := ExtractFQ([]string{"/nonexistent/path/reads.fq"}, false, "") if err == nil { t.Fatal("expected error for nonexistent input file, got nil") } } + +func writeTenXFastqExtract(t *testing.T, path string, ids []string, barcode string) { + t.Helper() + var sb strings.Builder + for _, id := range ids { + seq := barcode + "ACGTACGTAC" + sb.WriteString("@" + id + "\n" + seq + "\n+\n" + strings.Repeat("I", len(seq)) + "\n") + } + if err := os.WriteFile(path, []byte(sb.String()), 0o644); err != nil { + t.Fatal(err) + } +} + +func TestExtract10XFQ(t *testing.T) { + dir := t.TempDir() + r1 := filepath.Join(dir, "r1.fq") + wl := filepath.Join(dir, "whitelist.txt") + if err := os.WriteFile(wl, []byte("AAAACCCCGGGGTTTT\nTTTTGGGGCCCCAAAA\n"), 0o644); err != nil { + t.Fatal(err) + } + + var sb strings.Builder + sb.WriteString("@read1\nAAAACCCCGGGGTTTTACGTACGTAC\n+\nIIIIIIIIIIIIIIIIIIIIIIIIII\n") + sb.WriteString("@read2\nAAAACCCCGGGGTTTTACGTACGTAC\n+\nIIIIIIIIIIIIIIIIIIIIIIIIII\n") // duplicate barcode + sb.WriteString("@read3\nTTTTGGGGCCCCAAAAACGTACGTAC\n+\nIIIIIIIIIIIIIIIIIIIIIIIIII\n") // second barcode + sb.WriteString("@read4\nZZZZZZZZZZZZZZZZACGTACGTAC\n+\nIIIIIIIIIIIIIIIIIIIIIIIIII\n") // unrecognized -> dropped + if err := os.WriteFile(r1, []byte(sb.String()), 0o644); err != nil { + t.Fatal(err) + } + + var err error + out := captureStdout(t, func() { + err = ExtractFQ([]string{r1}, false, wl) + }) + if err != nil { + t.Fatalf("ExtractFQ (10X): %v", err) + } + lines := splitNonEmpty(string(out)) + if len(lines) != 2 { + t.Fatalf("got %d barcodes, want 2 (deduped, unrecognized dropped): %v", len(lines), lines) + } +} + +// ExtractFQ only ever reads fqs[0] for 10X input (the barcode lives solely +// in R1); passing a second file must not error, and its contents must be +// ignored rather than misread as a second R1. +func TestExtract10XFQ_IgnoresR2(t *testing.T) { + dir := t.TempDir() + r1 := filepath.Join(dir, "r1.fq") + r2 := filepath.Join(dir, "r2.fq") + wl := filepath.Join(dir, "whitelist.txt") + if err := os.WriteFile(wl, []byte("AAAACCCCGGGGTTTT\n"), 0o644); err != nil { + t.Fatal(err) + } + writeTenXFastqExtract(t, r1, []string{"read1"}, "AAAACCCCGGGGTTTT") + if err := os.WriteFile(r2, []byte("@read1\nGGGGCCCCAAAATTTT\n+\nIIIIIIIIIIIIIIII\n"), 0o644); err != nil { + t.Fatal(err) + } + + var err error + out := captureStdout(t, func() { + err = ExtractFQ([]string{r1, r2}, false, wl) + }) + if err != nil { + t.Fatalf("ExtractFQ (10X, paired): %v", err) + } + lines := splitNonEmpty(string(out)) + if len(lines) != 1 || lines[0] != "AAAACCCCGGGGTTTT" { + t.Errorf("got %v, want exactly [AAAACCCCGGGGTTTT]", lines) + } +} diff --git a/invalid/invalid_fq_test.go b/invalid/invalid_fq_test.go index cce28de..6d77fd5 100644 --- a/invalid/invalid_fq_test.go +++ b/invalid/invalid_fq_test.go @@ -1,8 +1,10 @@ package invalid import ( + "os" "path/filepath" "slices" + "strings" "testing" ) @@ -21,7 +23,7 @@ func TestFilterInvalidFQ_SplitsValidAndInvalid(t *testing.T) { prefix := filepath.Join(dir, "out") invPrefix := true - if err := FilterInvalidFQ([]string{r1}, prefix, invPrefix); err != nil { + if err := FilterInvalidFQ([]string{r1}, prefix, invPrefix, ""); err != nil { t.Fatalf("FilterInvalidFQ returned error: %v", err) } @@ -54,7 +56,7 @@ func TestFilterInvalidFQ_NoInvalidPrefix(t *testing.T) { writeHapFastq(t, r1, reads) prefix := filepath.Join(dir, "out") - if err := FilterInvalidFQ([]string{r1}, prefix, false); err != nil { + if err := FilterInvalidFQ([]string{r1}, prefix, false, ""); err != nil { t.Fatalf("FilterInvalidFQ returned error: %v", err) } @@ -78,7 +80,7 @@ func TestFilterInvalidFQ_PairedFiles(t *testing.T) { writeHapFastq(t, r2, []hapRead{{id: "p1", bc: "TTTT", vx: 1}}) prefix := filepath.Join(dir, "out") - if err := FilterInvalidFQ([]string{r1, r2}, prefix, false); err != nil { + if err := FilterInvalidFQ([]string{r1, r2}, prefix, false, ""); err != nil { t.Fatalf("FilterInvalidFQ returned error: %v", err) } if ids := readGzFastqIDs(t, prefix+".R1.fq.gz"); !slices.Contains(ids, "p1") { @@ -91,8 +93,137 @@ func TestFilterInvalidFQ_PairedFiles(t *testing.T) { func TestFilterInvalidFQ_MissingFile(t *testing.T) { dir := t.TempDir() - err := FilterInvalidFQ([]string{"/nonexistent/file.fq"}, filepath.Join(dir, "out"), false) + err := FilterInvalidFQ([]string{"/nonexistent/file.fq"}, filepath.Join(dir, "out"), false, "") if err == nil { t.Fatal("expected an error for a nonexistent input file") } } + +func writeTenXR1(t *testing.T, path string, reads [][2]string) { + t.Helper() + var sb strings.Builder + for _, r := range reads { + id, barcode := r[0], r[1] + seq := barcode + "ACGTACGTAC" + sb.WriteString("@" + id + "\n" + seq + "\n+\n" + strings.Repeat("I", len(seq)) + "\n") + } + if err := os.WriteFile(path, []byte(sb.String()), 0o644); err != nil { + t.Fatal(err) + } +} + +func writeTenXR2(t *testing.T, path string, ids []string) { + t.Helper() + var sb strings.Builder + for _, id := range ids { + sb.WriteString("@" + id + "\nGGGGCCCCAAAATTTT\n+\nIIIIIIIIIIIIIIII\n") + } + if err := os.WriteFile(path, []byte(sb.String()), 0o644); err != nil { + t.Fatal(err) + } +} + +// Regression test: FilterInvalid10XFQ's paired branch used to never write +// a valid R2 record to its own output file at all (only the keepInvalid +// path wrote R2, for invalid records) -- prefix.R2.fq.gz came out empty +// for every properly-paired, valid-barcode input. +func TestFilterInvalid10XFQ_PairedWritesValidR2(t *testing.T) { + dir := t.TempDir() + r1 := filepath.Join(dir, "r1.fq") + r2 := filepath.Join(dir, "r2.fq") + wl := filepath.Join(dir, "whitelist.txt") + if err := os.WriteFile(wl, []byte("AAAACCCCGGGGTTTT\n"), 0o644); err != nil { + t.Fatal(err) + } + writeTenXR1(t, r1, [][2]string{{"good1", "AAAACCCCGGGGTTTT"}}) + writeTenXR2(t, r2, []string{"good1"}) + + prefix := filepath.Join(dir, "out") + if err := FilterInvalidFQ([]string{r1, r2}, prefix, false, wl); err != nil { + t.Fatalf("FilterInvalidFQ (10X paired): %v", err) + } + + r1IDs := readGzFastqIDs(t, prefix+".R1.fq.gz") + r2IDs := readGzFastqIDs(t, prefix+".R2.fq.gz") + if !slices.Contains(r1IDs, "good1:AAAACCCCGGGGTTTT") { + t.Errorf("expected good1 in R1 output, got %v", r1IDs) + } + if len(r2IDs) == 0 { + t.Fatal("R2 output is empty -- valid R2 records must be written, not just invalid ones") + } + if !slices.Contains(r2IDs, "good1:AAAACCCCGGGGTTTT") { + t.Errorf("expected good1 in R2 output, got %v", r2IDs) + } +} + +// Regression test: FilterInvalid10XFQ used to have no separate invalid-output +// file at all for 10X input, unlike every other format -- "kept" invalid +// records were written into the same file as valid ones instead of their +// own prefix.invalid.R*.fq.gz. +func TestFilterInvalid10XFQ_SeparatesInvalidOutput(t *testing.T) { + dir := t.TempDir() + r1 := filepath.Join(dir, "r1.fq") + r2 := filepath.Join(dir, "r2.fq") + wl := filepath.Join(dir, "whitelist.txt") + if err := os.WriteFile(wl, []byte("AAAACCCCGGGGTTTT\n"), 0o644); err != nil { + t.Fatal(err) + } + writeTenXR1(t, r1, [][2]string{ + {"good1", "AAAACCCCGGGGTTTT"}, + {"bad1", "ZZZZZZZZZZZZZZZZ"}, + }) + writeTenXR2(t, r2, []string{"good1", "bad1"}) + + prefix := filepath.Join(dir, "out") + if err := FilterInvalidFQ([]string{r1, r2}, prefix, true, wl); err != nil { + t.Fatalf("FilterInvalidFQ (10X paired, keepInvalid): %v", err) + } + + validR1 := readGzFastqIDs(t, prefix+".R1.fq.gz") + validR2 := readGzFastqIDs(t, prefix+".R2.fq.gz") + invalidR1 := readGzFastqIDs(t, prefix+".invalid.R1.fq.gz") + invalidR2 := readGzFastqIDs(t, prefix+".invalid.R2.fq.gz") + + for _, ids := range [][]string{validR1, validR2} { + if slices.ContainsFunc(ids, func(s string) bool { return strings.HasPrefix(s, "bad1") }) { + t.Errorf("invalid record bad1 leaked into the main valid output: %v", ids) + } + } + for _, ids := range [][]string{invalidR1, invalidR2} { + if !slices.ContainsFunc(ids, func(s string) bool { return strings.HasPrefix(s, "bad1") }) { + t.Errorf("expected bad1 in the separate invalid output, got %v", ids) + } + if slices.ContainsFunc(ids, func(s string) bool { return strings.HasPrefix(s, "good1") }) { + t.Errorf("valid record good1 leaked into the invalid output: %v", ids) + } + } +} + +func TestFilterInvalid10XFQ_SingleEnd(t *testing.T) { + dir := t.TempDir() + r1 := filepath.Join(dir, "r1.fq") + wl := filepath.Join(dir, "whitelist.txt") + if err := os.WriteFile(wl, []byte("AAAACCCCGGGGTTTT\n"), 0o644); err != nil { + t.Fatal(err) + } + writeTenXR1(t, r1, [][2]string{ + {"good1", "AAAACCCCGGGGTTTT"}, + {"bad1", "ZZZZZZZZZZZZZZZZ"}, + }) + + prefix := filepath.Join(dir, "out") + if err := FilterInvalidFQ([]string{r1}, prefix, true, wl); err != nil { + t.Fatalf("FilterInvalidFQ (10X single-end, keepInvalid): %v", err) + } + if _, err := os.Stat(prefix + ".R2.fq.gz"); err == nil { + t.Error("single-end 10X input must not produce an R2 output file") + } + validIDs := readGzFastqIDs(t, prefix+".R1.fq.gz") + invalidIDs := readGzFastqIDs(t, prefix+".invalid.R1.fq.gz") + if !slices.ContainsFunc(validIDs, func(s string) bool { return strings.HasPrefix(s, "good1") }) { + t.Errorf("expected good1 in valid output, got %v", validIDs) + } + if !slices.ContainsFunc(invalidIDs, func(s string) bool { return strings.HasPrefix(s, "bad1") }) { + t.Errorf("expected bad1 in invalid output, got %v", invalidIDs) + } +} diff --git a/stats/stats_fq_test.go b/stats/stats_fq_test.go index fade62b..60b0134 100644 --- a/stats/stats_fq_test.go +++ b/stats/stats_fq_test.go @@ -24,7 +24,7 @@ func writeFastq(t *testing.T, records ...string) string { func statsFQOf(t *testing.T, path, platform string) string { t.Helper() - out, err := captureStdout(t, func() error { return StatsFQ(path, platform) }) + out, err := captureStdout(t, func() error { return StatsFQ(path, platform, "") }) if err != nil { t.Fatalf("StatsFQ: %v", err) } @@ -156,7 +156,7 @@ func TestStatsFQ_StlfrAndTellseqByName(t *testing.T) { func TestStatsFQ_Errors(t *testing.T) { noBarcodes := writeFastq(t, fq("r1 1:N:0:ATGGAGA"), fq("r2 1:N:0:ATGGAGA")) - if _, err := captureStdout(t, func() error { return StatsFQ(noBarcodes, "") }); err == nil || + if _, err := captureStdout(t, func() error { return StatsFQ(noBarcodes, "", "") }); err == nil || !strings.Contains(err.Error(), "unable to determine") { t.Errorf("no detectable technology: err = %v", err) } @@ -166,11 +166,45 @@ func TestStatsFQ_Errors(t *testing.T) { t.Errorf("no barcodes with a named platform: got:\n%s\nwant:\n%s", got, want) } - if _, err := captureStdout(t, func() error { return StatsFQ(noBarcodes, "10x") }); err == nil || + if _, err := captureStdout(t, func() error { return StatsFQ(noBarcodes, "10x", "") }); err == nil || !strings.Contains(err.Error(), "unknown platform") { t.Errorf("unknown platform: err = %v", err) } - if err := StatsFQ(fmt.Sprintf("%s/missing.fq", t.TempDir()), Haplotagging); err == nil { + if err := StatsFQ(fmt.Sprintf("%s/missing.fq", t.TempDir()), Haplotagging, ""); err == nil { t.Error("a missing file should be an error") } } + +func tenxRecord(id, barcode string) string { + seq := barcode + "ACGTACGTAC" + return "@" + id + "\n" + seq + "\n+\n" + strings.Repeat("I", len(seq)) + "\n" +} + +// Regression test: StatsFQ's 10X branch built parse/sch from bcMap but then +// unconditionally called chooseParser(fq, platform) right after, which +// reassigned parse/sch every time regardless of bcMap -- --bc-list was a +// complete no-op and 10X input always fell through to plain +// haplotagging/stlfr/tellseq auto-detection (which would just fail, since +// 10X has no self-describing marker). +func TestStatsFQ_TenX(t *testing.T) { + dir := t.TempDir() + wl := fmt.Sprintf("%s/whitelist.txt", dir) + if err := os.WriteFile(wl, []byte("AAAACCCCGGGGTTTT\n"), 0o644); err != nil { + t.Fatal(err) + } + + path := writeFastq(t, + tenxRecord("read1", "AAAACCCCGGGGTTTT"), + tenxRecord("read2", "AAAACCCCGGGGTTTT"), + tenxRecord("read3", "ZZZZZZZZZZZZZZZZ"), // not in whitelist -> no barcode at all + ) + + got, err := captureStdout(t, func() error { return StatsFQ(path, "", wl) }) + if err != nil { + t.Fatalf("StatsFQ (10X): %v", err) + } + lines := strings.Split(strings.TrimRight(got, "\n"), "\n") + if lines[0] != "Reads\t3" || lines[1] != "Barcodes\t2" || lines[2] != "Barcodes_Valid\t2" || lines[3] != "Barcodes_Invalid\t0" { + t.Errorf("unexpected 10X stats output:\n%s", got) + } +} diff --git a/unlinked/unlinked_fq_test.go b/unlinked/unlinked_fq_test.go index e1b7dbc..076ceb8 100644 --- a/unlinked/unlinked_fq_test.go +++ b/unlinked/unlinked_fq_test.go @@ -4,6 +4,7 @@ import ( "os" "path/filepath" "slices" + "strings" "testing" ) @@ -88,9 +89,8 @@ func TestRemoveUnlinkedFQ_SplitsSingletonsAndNot(t *testing.T) { writeHapFastq(t, r1, reads) prefix := filepath.Join(dir, "out") - bcCountFile := filepath.Join(dir, "counts.tsv") - if err := RemoveUnlinkedFQ([]string{r1}, prefix, true, bcCountFile); err != nil { + if err := RemoveUnlinkedFQ([]string{r1}, prefix, true, true, ""); err != nil { t.Fatalf("RemoveUnlinkedFQ returned error: %v", err) } @@ -110,7 +110,7 @@ func TestRemoveUnlinkedFQ_SplitsSingletonsAndNot(t *testing.T) { t.Errorf("invalid-barcode read g1 should be dropped entirely, got main=%v single=%v", mainIDs, singleIDs) } - if _, err := os.Stat(bcCountFile); err != nil { + if _, err := os.Stat(prefix + ".count"); err != nil { t.Errorf("expected barcode count file to be created: %v", err) } } @@ -126,7 +126,7 @@ func TestRemoveUnlinkedFQ_NoSingletonOutput(t *testing.T) { writeHapFastq(t, r1, reads) prefix := filepath.Join(dir, "out") - if err := RemoveUnlinkedFQ([]string{r1}, prefix, false, ""); err != nil { + if err := RemoveUnlinkedFQ([]string{r1}, prefix, false, false, ""); err != nil { t.Fatalf("RemoveUnlinkedFQ returned error: %v", err) } mainIDs := readGzFastqIDs(t, prefix+".R1.fq.gz") @@ -137,8 +137,109 @@ func TestRemoveUnlinkedFQ_NoSingletonOutput(t *testing.T) { func TestRemoveUnlinkedFQ_MissingFile(t *testing.T) { dir := t.TempDir() - err := RemoveUnlinkedFQ([]string{"/nonexistent/file.fq"}, filepath.Join(dir, "out"), false, "") + err := RemoveUnlinkedFQ([]string{"/nonexistent/file.fq"}, filepath.Join(dir, "out"), false, false, "") if err == nil { t.Fatal("expected an error for a nonexistent input file") } } + +func writeTenXR1Unlinked(t *testing.T, path string, reads [][2]string) { + t.Helper() + var sb strings.Builder + for _, r := range reads { + id, barcode := r[0], r[1] + seq := barcode + "ACGTACGTAC" + sb.WriteString("@" + id + "\n" + seq + "\n+\n" + strings.Repeat("I", len(seq)) + "\n") + } + if err := os.WriteFile(path, []byte(sb.String()), 0o644); err != nil { + t.Fatal(err) + } +} + +func writeTenXR2Unlinked(t *testing.T, path string, ids []string) { + t.Helper() + var sb strings.Builder + for _, id := range ids { + sb.WriteString("@" + id + "\nGGGGCCCCAAAATTTT\n+\nIIIIIIIIIIIIIIII\n") + } + if err := os.WriteFile(path, []byte(sb.String()), 0o644); err != nil { + t.Fatal(err) + } +} + +// Regression test: RemoveUnlinkedFQ accepted a bcmap parameter but never +// used it anywhere in the function body, so --bc-list/-b was silently a +// no-op for rm-unlinked -- a genuine 10X input would always fall through +// to auto-detection and fail with "unable to determine linked-read +// technology", since 10X has no self-describing marker to detect. +func TestRemoveUnlinked10XFQ_SplitsSingletonsAndNot(t *testing.T) { + dir := t.TempDir() + r1 := filepath.Join(dir, "r1.fq") + r2 := filepath.Join(dir, "r2.fq") + wl := filepath.Join(dir, "whitelist.txt") + if err := os.WriteFile(wl, []byte("AAAACCCCGGGGTTTT\nTTTTGGGGCCCCAAAA\n"), 0o644); err != nil { + t.Fatal(err) + } + + writeTenXR1Unlinked(t, r1, [][2]string{ + {"a1", "AAAACCCCGGGGTTTT"}, + {"a2", "AAAACCCCGGGGTTTT"}, // same barcode as a1 -> linked + {"c1", "TTTTGGGGCCCCAAAA"}, // singleton + }) + writeTenXR2Unlinked(t, r2, []string{"a1", "a2", "c1"}) + + prefix := filepath.Join(dir, "out") + if err := RemoveUnlinkedFQ([]string{r1, r2}, prefix, true, true, wl); err != nil { + t.Fatalf("RemoveUnlinkedFQ (10X): %v", err) + } + + mainR1 := readGzFastqIDs(t, prefix+".R1.fq.gz") + mainR2 := readGzFastqIDs(t, prefix+".R2.fq.gz") + singleR1 := readGzFastqIDs(t, prefix+".unlinked.R1.fq.gz") + singleR2 := readGzFastqIDs(t, prefix+".unlinked.R2.fq.gz") + + hasPrefix := func(ids []string, id string) bool { + return slices.ContainsFunc(ids, func(s string) bool { return strings.HasPrefix(s, id) }) + } + + if !hasPrefix(mainR1, "a1") || !hasPrefix(mainR1, "a2") { + t.Errorf("expected a1,a2 in main R1 output, got %v", mainR1) + } + if !hasPrefix(mainR2, "a1") || !hasPrefix(mainR2, "a2") { + t.Errorf("expected a1,a2 in main R2 output, got %v", mainR2) + } + if hasPrefix(mainR1, "c1") || hasPrefix(mainR2, "c1") { + t.Errorf("singleton c1 should not be in main output, got R1=%v R2=%v", mainR1, mainR2) + } + if !hasPrefix(singleR1, "c1") || !hasPrefix(singleR2, "c1") { + t.Errorf("expected c1 in singleton output, got R1=%v R2=%v", singleR1, singleR2) + } + if _, err := os.Stat(prefix + ".count"); err != nil { + t.Errorf("expected barcode count file to be created: %v", err) + } +} + +func TestRemoveUnlinked10XFQ_SingleEnd(t *testing.T) { + dir := t.TempDir() + r1 := filepath.Join(dir, "r1.fq") + wl := filepath.Join(dir, "whitelist.txt") + if err := os.WriteFile(wl, []byte("AAAACCCCGGGGTTTT\n"), 0o644); err != nil { + t.Fatal(err) + } + writeTenXR1Unlinked(t, r1, [][2]string{ + {"a1", "AAAACCCCGGGGTTTT"}, + {"a2", "AAAACCCCGGGGTTTT"}, + }) + + prefix := filepath.Join(dir, "out") + if err := RemoveUnlinkedFQ([]string{r1}, prefix, false, false, wl); err != nil { + t.Fatalf("RemoveUnlinkedFQ (10X single-end): %v", err) + } + if _, err := os.Stat(prefix + ".R2.fq.gz"); err == nil { + t.Error("single-end 10X input must not produce an R2 output file") + } + mainR1 := readGzFastqIDs(t, prefix+".R1.fq.gz") + if len(mainR1) != 2 { + t.Errorf("got %d records in main output, want 2: %v", len(mainR1), mainR1) + } +} From 7067927fe018fe530ba50b24de27912292e496b2 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 21 Oct 2026 18:35:35 +0000 Subject: [PATCH 5/6] docs: fix stale --bc-map references, document new 10X support convert.mdx, sort.mdx, and ncbi.mdx still named the flag --bc-map, which was renamed to --bc-list; fixed throughout, and documented the "10x"/"10X" embedded-whitelist sentinel (fastq.Load10X) that none of the docs mentioned yet. count.mdx, extract.mdx, rm-invalid.mdx, and rm-unlinked.mdx gained real 10X support in the same commit as this docs update but had no mention of it at all -- added a --bc-list row to each options table and a short paragraph explaining why 10X needs it (no self-describing marker, so it can't be auto-detected like haplotagging/stlfr/tellseq) and, for count/extract, that only the first (R1) file is read. rm-unlinked.mdx's options table also still described --bc-count as a string flag with the (now reassigned) -b shorthand; it's a bool that writes to PREFIX.count, documented accordingly. Note: `djinn fq stats` still has no docs page at all (10X-related or otherwise) -- out of scope here since it's not a 10X-specific gap, but worth a follow-up. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_015tF7ppXiNo2Adbv8YXcwGL --- docs/src/content/docs/fastq/convert.mdx | 15 +++++++++------ docs/src/content/docs/fastq/count.mdx | 7 +++++++ docs/src/content/docs/fastq/extract.mdx | 7 +++++++ docs/src/content/docs/fastq/ncbi.mdx | 8 +++++++- docs/src/content/docs/fastq/rm-invalid.mdx | 8 ++++++++ docs/src/content/docs/fastq/rm-unlinked.mdx | 8 +++++++- docs/src/content/docs/fastq/sort.mdx | 7 +++++-- 7 files changed, 50 insertions(+), 10 deletions(-) diff --git a/docs/src/content/docs/fastq/convert.mdx b/docs/src/content/docs/fastq/convert.mdx index 56b8753..d6f5334 100644 --- a/docs/src/content/docs/fastq/convert.mdx +++ b/docs/src/content/docs/fastq/convert.mdx @@ -9,10 +9,13 @@ import { Badge } from '@astrojs/starlight/components'; If you need to convert between the linked-read technologies' distinct FASTQ formats and their -barcoding systems, use `convert`. It auto-detects the input format as one of 10X, haplotagging, TELLseq, or stLFR, -and converts it to the format provided as the `TARGET` positional argument. If the -input data is in 10X format, where the barcode is usually the first 16 bases in the R1 sequence, -you will need to provide a `--bc-map` file that lets Djinn know what barcodes to look for. +barcoding systems, use `convert`. It auto-detects the input format as one of haplotagging, TELLseq, or stLFR +by inspecting the first 100 records. 10X cannot be auto-detected this way (its barcode is just the +first 16 bases of the R1 sequence, with nothing marking it as a barcode), so converting **from** 10X +requires `--bc-list` to tell Djinn which 16-base sequences are real barcodes: either a path to your +own whitelist file, or the literal value `10x` to use the official 10X Genomics barcode whitelist +embedded in Djinn. Converting **to** 10X draws new barcodes from that same embedded whitelist and +needs no `--bc-list` at all. In all cases, a file will be created with the barcode conversion map. ## Usage @@ -36,11 +39,11 @@ convert 10x gcurratum_as10X gcurratum.R1.fq.gz gcurratum.R2.fq.gz | `PREFIX` | output filename prefix | | `TARGET` | output barcode and FASTQ style: `10x`, `haplotagging`, `tellseq`, `stlfr`, `standard` | | `file.R*.fq.gz` | FASTQ file or file pair | -| `-b` `--bc-map` | Valid barcode list | +| `-b` `--bc-list` | Valid barcode list (or literal `10x` for the embedded whitelist) |
- example `--bc-map` file + example `--bc-list` file ``` ATGGAAGCCGTAGTTA ACGGAAGCCGTAGTTC diff --git a/docs/src/content/docs/fastq/count.mdx b/docs/src/content/docs/fastq/count.mdx index 9796e53..93b3db6 100644 --- a/docs/src/content/docs/fastq/count.mdx +++ b/docs/src/content/docs/fastq/count.mdx @@ -13,6 +13,12 @@ You can use one file of a pair (e.g. R1 or R2), but including both files will pr and only count occurences of _new_ barcodes in the second file that aren't present in the first file. This is done to avoid double-counting barcodes. +Format is auto-detected from the first 100 records for haplotagging, TELLseq, and stLFR. 10X data +cannot be auto-detected this way (its barcode is just the first 16 bases of R1, with nothing marking +it as a barcode), so counting 10X barcodes requires `--bc-list`: either a path to your own +whitelist, or the literal value `10x` to use the official 10X Genomics whitelist embedded in Djinn. +For 10X input, only the first file (R1) is read, since that's the only place the barcode lives. + ## Usage ```bash @@ -29,3 +35,4 @@ djinn fastq count sample1.F.fq.gz sample1.R.fq.gz > sample1.bxcount |:------------------|:------------------------------------------------------------------------------| | `file.R*.fq.gz` | FASTQ file or file pair | | `-i` `--invalid` | Include invalid barcodes | +| `-b` `--bc-list` | Valid barcode list, or literal `10x` for the embedded official whitelist | diff --git a/docs/src/content/docs/fastq/extract.mdx b/docs/src/content/docs/fastq/extract.mdx index 01b40d4..782a91e 100644 --- a/docs/src/content/docs/fastq/extract.mdx +++ b/docs/src/content/docs/fastq/extract.mdx @@ -8,6 +8,12 @@ import { Badge } from '@astrojs/starlight/components'; To get a list of all the unique barcodes in your data, use `extract`. Writes to `stdout`. +Format is auto-detected from the first 100 records for haplotagging, TELLseq, and stLFR. 10X data +cannot be auto-detected this way (its barcode is just the first 16 bases of R1, with nothing marking +it as a barcode), so extracting 10X barcodes requires `--bc-list`: either a path to your own +whitelist, or the literal value `10x` to use the official 10X Genomics whitelist embedded in Djinn. +For 10X input, only the first file (R1) is read, since that's the only place the barcode lives. + ## Usage ```bash djinn fastq extract file.R1.fq.gz [file.R2.fq.gz] @@ -23,3 +29,4 @@ djinn fastq extract -i sample1.F.fq.gz sample1.R.fq.gz > sample1.bc |:------------------|:------------------------------------------------------------------------------| | `file.R*fq.gz` | FASTQ file or file pair | | `-i` `--invalid` | Include invalid barcodes | +| `-b` `--bc-list` | Valid barcode list, or literal `10x` for the embedded official whitelist | diff --git a/docs/src/content/docs/fastq/ncbi.mdx b/docs/src/content/docs/fastq/ncbi.mdx index 2b69309..2a2c4f2 100644 --- a/docs/src/content/docs/fastq/ncbi.mdx +++ b/docs/src/content/docs/fastq/ncbi.mdx @@ -23,6 +23,12 @@ files into a single unaligned BAM file without any data loss (and also free up d you risk losing the barcode forever! +Format is auto-detected from the first 100 records for haplotagging, TELLseq, and stLFR. 10X cannot be +auto-detected this way (its barcode is just the first 16 bases of R1, with nothing marking it as a +barcode), so converting 10X FASTQ to BAM requires `--bc-list` to tell Djinn which 16-base sequences +are real barcodes -- either a path to your own whitelist, or the literal value `10x` to use the +official 10X Genomics whitelist embedded in Djinn. + ## Usage ```bash @@ -43,5 +49,5 @@ samtools import -O BAM -T "*" -1 file.R1.fq file.R2.fq > out.bam | argument | description | |:-----------------|:---------------------------------------------------------------------------| | `file.R*.fq.gz` | FASTQ file or file pair | -| `-b` `--bc-map` | Valid barcode list | +| `-b` `--bc-list` | Valid barcode list, or literal `10x` for the embedded official whitelist | | `-t` `--threads` | Threads to use (default: 2) | diff --git a/docs/src/content/docs/fastq/rm-invalid.mdx b/docs/src/content/docs/fastq/rm-invalid.mdx index 832dc81..f57b03c 100644 --- a/docs/src/content/docs/fastq/rm-invalid.mdx +++ b/docs/src/content/docs/fastq/rm-invalid.mdx @@ -10,6 +10,13 @@ You may be interested in retaining only "useful" linked reads, that is, reads wi Using `rm-invalid` will remove alignment records invalid barcodes. If you want to keep a separate file (or files) of records with invalid barcodes use the `--invalid` flag. +Format is auto-detected from the first 100 records for haplotagging, TELLseq, and stLFR. 10X data +cannot be auto-detected this way (its barcode is just the first 16 bases of R1, with nothing marking +it as a barcode), so filtering 10X data requires `--bc-list`: either a path to your own whitelist, +or the literal value `10x` to use the official 10X Genomics whitelist embedded in Djinn. A 10X R1 +read whose first 16 bases don't match any whitelist entry has no barcode at all and is treated as +invalid. + ## Usage ```bash djinn sam rm-invalid [--invalid] PREFIX file.R1.fq.gz [file.R2.fq.gz] @@ -26,3 +33,4 @@ djinn sam rm-invalid --invalid echinda.valid echidna.R*.fq.gz | `PREFIX` | output filename prefix | | `file.R*.fq.gz` | FASTQ file or file pair | | `-i` `--invalid` | output records with invalid barcodes to `PREFIX.invalid.*` | +| `-b` `--bc-list` | Valid barcode list, or literal `10x` for the embedded official whitelist | diff --git a/docs/src/content/docs/fastq/rm-unlinked.mdx b/docs/src/content/docs/fastq/rm-unlinked.mdx index e5d770a..2de7e5e 100644 --- a/docs/src/content/docs/fastq/rm-unlinked.mdx +++ b/docs/src/content/docs/fastq/rm-unlinked.mdx @@ -12,6 +12,11 @@ to isolate records with meaningful linked-read information. If you want to keep file (or files) of unlinked/singleton records, use the `--unlinked` flag. Barcodes are counted per read-pair or unpair read. +Format is auto-detected from the first 100 records for haplotagging, TELLseq, and stLFR. 10X data +cannot be auto-detected this way (its barcode is just the first 16 bases of R1, with nothing marking +it as a barcode), so filtering 10X data requires `--bc-list`: either a path to your own whitelist, +or the literal value `10x` to use the official 10X Genomics whitelist embedded in Djinn. + ## Usage ```bash @@ -28,5 +33,6 @@ djinn fastq rm-unlinked --unlinked hyena.clean hyena.R*.fq.gz | :------------------ | :--------------------------------------------------------------------------------- | | `PREFIX` | output filename prefix | | `file.R*.fq.gz` | FASTQ file or file pair | -| `-b` `--bc-count` | write valid barcode counts to this file | +| `--bc-count` | write valid barcode counts to `PREFIX.count` | | `-u` `--unlinked` | write valid singleton records to `PREFIX.unlinked` | +| `-b` `--bc-list` | Valid barcode list, or literal `10x` for the embedded official whitelist | diff --git a/docs/src/content/docs/fastq/sort.mdx b/docs/src/content/docs/fastq/sort.mdx index 37c6dac..fff11f0 100644 --- a/docs/src/content/docs/fastq/sort.mdx +++ b/docs/src/content/docs/fastq/sort.mdx @@ -14,7 +14,10 @@ in the process. Sorting is done using a `samtools` style k-way merge, where reads (or read pairs) are inlined with an intermediate format and cast across several temporary files located at `--tmp-dir`, where -they are then merged back together. Using the `--standardize` flag will result in sorted +they are then merged back together. Format is auto-detected from the first 100 records for +haplotagging, TELLseq, and stLFR; 10X cannot be auto-detected this way (its barcode is just the +first 16 bases of R1, with nothing marking it as a barcode), so sorting 10X data requires +`--bc-list` to tell Djinn which 16-base sequences are real barcodes. Using the `--standardize` flag will result in sorted [standard format](https://blinkseq.github.io/lastq/) data that, among other things, is ready for use in [arachne](https://github.com/blinkseq/arachne). @@ -40,7 +43,7 @@ When input is a FASTQ file pair, the reads are expected to be synced between R1 |:------------------|:------------------------------------------------------------------------------| | `PREFIX` | output filename prefix | | `file.R*.fq.gz` | FASTQ file or file pair | -| `-b` `--bc-map` | valid barcode list | +| `-b` `--bc-list` | valid barcode list, or literal `10x` for the embedded official whitelist | | `-s` `--standardize` | output using standardized BX:Z:/VX:i: tags instead of the detected input format | | `-@` `--threads` | treads to use (default: 2) | | `-t` `--tmp-prefix` | folder for temporary files | From 47a5b20296241f6d9f222428f8e703023098fff0 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 21 Oct 2026 18:41:03 +0000 Subject: [PATCH 6/6] docs: add missing stats pages for fq and sam Both `djinn fq stats` and `djinn sam stats` had no documentation page at all. Added one for each, matching the existing pages' format, tone, and length: a short intro, a Usage section with an example, an Options table, and (new here, since neither existing page needed it) an Output table describing the columns each command actually prints. fastq/stats.mdx also documents the --bc-list flag and the 10X embedded-whitelist sentinel, consistent with the other FASTQ docs pages updated for 10X support. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_015tF7ppXiNo2Adbv8YXcwGL --- docs/src/content/docs/fastq/stats.mdx | 49 ++++++++++++++++++++++++ docs/src/content/docs/sam/stats.mdx | 54 +++++++++++++++++++++++++++ 2 files changed, 103 insertions(+) create mode 100644 docs/src/content/docs/fastq/stats.mdx create mode 100644 docs/src/content/docs/sam/stats.mdx diff --git a/docs/src/content/docs/fastq/stats.mdx b/docs/src/content/docs/fastq/stats.mdx new file mode 100644 index 0000000..95ea170 --- /dev/null +++ b/docs/src/content/docs/fastq/stats.mdx @@ -0,0 +1,49 @@ +--- +title: Barcode statistics for a FASTQ file +sidebar: + order: 4 + label: stats +--- +import { Badge } from '@astrojs/starlight/components'; + +If you want a quick read on how clean your linked-read barcodes are, use `stats`. It reports how many +reads carry a barcode, how many of those barcodes are valid, and, for the barcodes that aren't, which +position tends to go wrong. Writes to `stdout`. + +Format is auto-detected from the first 100 records for haplotagging, TELLseq, and stLFR, or you can +name it with `--platform` to skip detection. 10X data cannot be auto-detected this way (its barcode is +just the first 16 bases of R1, with nothing marking it as a barcode), so 10X input requires `--bc-list` +instead: either a path to your own whitelist, or the literal value `10x` to use the official 10X +Genomics whitelist embedded in Djinn. + +## Usage +```bash +djinn fastq stats [options] file.fq.gz +``` + +Example +```bash +djinn fastq stats curratum.R1.fq.gz > curratum.R1.stats +``` + +## Options +| argument | description | +|:--------------------|:----------------------------------------------------------------------------------------------| +| `file.fq.gz` | FASTQ file | +| `-p` `--platform` | Linked-read technology: `haplotagging`, `stlfr` or `tellseq` (default: detect from the first 100 reads) | +| `-b` `--bc-list` | Valid barcode list, or literal `10x` for the embedded official whitelist | +| `-@` `--threads` | Threads to use (default: 2) | + +## Output + +| column | description | +|:-------------------|:-------------------------------------------------------------------------| +| `Reads` | total sequences in the file | +| `Barcodes` | reads that carry a linked-read barcode | +| `Barcodes_Valid` | barcodes that are valid | +| `Barcodes_Invalid` | barcodes that are not valid | + +followed by one line per barcode position that can go wrong, with the number of barcodes that were bad +at that position: haplotagging segments `A`, `C`, `B`, `D` (e.g. a segment reading `00`), stLFR parts +`1`-`3` (a part reading `0`), TELLseq bases `1`-`18` (an `N`), or 10X bases `1`-`16` (an `N`). A barcode +with several bad positions adds to each of them. diff --git a/docs/src/content/docs/sam/stats.mdx b/docs/src/content/docs/sam/stats.mdx new file mode 100644 index 0000000..239371a --- /dev/null +++ b/docs/src/content/docs/sam/stats.mdx @@ -0,0 +1,54 @@ +--- +title: Linked-read metrics per molecule +sidebar: + order: 4 + label: stats +--- +import { Badge, Aside } from '@astrojs/starlight/components'; + +Once your reads are aligned, `stats` reports per-molecule metrics: how many reads a barcode has on a +contig, where they start and end, how many bases aligned, and how well the molecule is covered. A +"molecule" is one barcode's reads on one contig; use `--distance-threshold` to split a barcode's reads +into separate molecules (`BARCODE-1`, `BARCODE-2`, ...) when the gap between consecutive reads is too +large, the same splitting `deconvolve` does, without writing a new BAM. Writes a table to `stdout`. + +Barcodes come from the `BX:Z` tag (or a TELLseq/stLFR-style read name), and a barcode is valid if the +`VX:i` tag says so, or, when there is no `VX:i` tag, if the barcode itself is valid. Reads without a +valid barcode are pooled into one `invalid` line per contig. + + + +## Usage +```bash +djinn sam stats [options] file.bam +``` + +Example +```bash +djinn sam stats -d 50000 -@ 4 curratum.bam > curratum.molecules.tsv +``` + +## Options +| argument | description | +|:------------------------------|:------------------------------------------------------------------------------------------------------| +| `file.bam` | coordinate-sorted SAM/BAM file, or `-` for stdin | +| `-d` `--distance-threshold` | distance (bp) that splits a barcode's reads into separate molecules, where `0` skips splitting (default: 50000) | +| `-@` `--threads` | Threads to use (default: 2) | + +## Output + +| column | description | +|:-------------------|:------------------------------------------------------------------------| +| `contig` | chromosome/contig the molecule is on | +| `molecule` | the barcode, suffixed with `-N` if split by `--distance-threshold` | +| `reads` | number of reads sharing the barcode on this molecule | +| `start` | leftmost aligned position | +| `end` | rightmost aligned position | +| `length_inferred` | `end - start` | +| `aligned_bp` | total aligned bases across the molecule's reads | +| `insert_len` | inferred insert length | +| `coverage_bp` | `aligned_bp` / `length_inferred` | +| `coverage_inserts` | inferred-insert coverage of the molecule |