Repository navigation
First release of LRSinglecell. The pipeline is a continuation of the scdnalong branch of
ljwharbers/scnanoseq (at 1ff784a), itself built on
nf-core/scnanoseq v1.3.0, rebuilt on the nf-core/tools 4.2.1 template as a
standalone IntGenomicsLab pipeline. The changes below are relative to nf-core/scnanoseq v1.3.0.
Pipeline
- Renamed to
IntGenomicsLab/lrsinglecell; the entry workflow isINTGENOMICSLAB_LRSINGLECELLand the software versions file ispipeline_info/lrsinglecell_software_mqc_versions.yml - Test profiles read their data from IntGenomicsLab/test-datasets (
lrsinglecell/), copied unchanged from nf-core/test-datasetsscnanoseq - Moved to the nf-core/tools 4.2.1 template (minimum Nextflow
25.10.4,nf-schema2.7.3,nft-utils1.2.1, updatedutils_*subworkflows,multiqcmodule andfastqcmodule (FastQC 0.13.0, only used with--skip_fastqc false))
Enhancements
- DNA samples are now deduplicated by
mark_dna_duplicates.pyinstead ofpicard MarkDuplicates. Within a barcode it keeps Picard's single-end key and also merges the two records flexiplex writes for one raw read, up to 10 bp of 5' jitter between reads that read through to the far adapter, and A-A fragments read from the other end; onchrM/MTonly the exact relations apply. On three libraries this recovers 99% of the same-molecule pairs Picard missed (judged by heterozygous-SNP haplotypes) without adding false merges, raising the DNA duplicate rate by 1-2.5 points. The duplicate flag now also reaches secondary and supplementary alignments, and the marked BAM is published as written, without a second sort. Metrics stay in Picard format for MultiQC, next to a per-rule*.dedup_summary.tsv; seedocs/output.md - Added support for processing single-cell/nuclei DNA samples alongside cDNA samples in the same pipeline run, via a new optional
typecolumn (dna/cdna) in the input samplesheet - Added
flexiplex/discovery,flexiplex/filterandflexiplex/assignmodules for barcode extraction, filtering and assignment - Added
demultiplex_flexiplexanddemultiplex_blazesubworkflows, replacing the previous inline demultiplexing steps in the main workflow - Added
align_deduplicate_dnasubworkflow (minimap2alignment,mark_dna_duplicates.py,BAM_STATS_SAMTOOLS) for DNA sample processing - Added
demux_tool_cdna/demux_tool_dnaparameters to selectflexiplexorblazefor cDNA demultiplexing (DNA currently supportsflexiplexonly) - Added
lrkallistoas a third--quantifieroption: an alignment-free path that pseudoaligns the demultiplexed FASTQ with lr-kallisto (kallisto bus --long), counts genes withbustools count --genecountsand quantifies transcripts withkallisto quant-tcc --long, deduplicating UMIs withbustoolsin both cases. Adds atest_lrkallistoprofile and thegffreadmodule for transcript sequence extraction - Bumped IsoQuant to v3.13.1
- Deduplication now groups UMIs by feature rather than by alignment position. Genome alignments are tagged with a gene by gene-body overlap (new
tag_genesandsplit_gene_statusmodules,GX/GN/GStags) and grouped withumi_tools --per-gene; reads with an ambiguous gene call or none are deduplicated positionally and merged back, so no reads are dropped. Transcriptome alignments useumi_tools --per-gene --per-contig, where a contig is a transcript. This substantially changes absolute UMI counts — on the data this was developed against, total duplicate removal rose from 6.5% to ~34%, i.e. previous counts were inflated by roughly 1.6x. Relative expression is essentially unchanged (Spearman 0.998, no genes lost), so clustering and ratio-based analyses are unaffected, but UMIs per cell, library complexity and saturation estimates will differ from earlier releases. Revert with--dedup_per_gene=false/--dedup_per_contig=false - Added a per-sample UMI deduplication summary (
<sample>/<type>/{genome,transcriptome}/qc/umitools/dedup/*.umi_dedup_summary.tsv) with input reads, output reads, duplicate reads and duplication rate, summed over the per-chunkumi_tools dedupruns by the newumitools_dedup_summarymodule, and shown as a "UMI-tools Deduplication" table in the MultiQC report. Previously theumi_toolslogs were discarded and no duplication rate was reported umi_tools dedupnow always takes the cell barcode from the correctedCBtag. Previously this applied only when--demux_tool_cdna flexiplexwas set; theblazepath fell back to reading the raw, uncorrected barcode out of the read name, which fragmented cells and made the two demultiplexing paths incomparable- Bumped flexiplex to v1.02.7 and run it with
-a true, so it reports every read rather than only the reads it could assign. Reads from a droplet below the knee are no longer dropped for want of a known barcode, and the uncorrected barcode comes back inCR, which is what makes such a droplet distinguishable from a read with no barcode region at all. v1.02.7 also fixes an out-of-bounds read in flexiplex's UMI handling that only triggers on UMI-less protocols, i.e. the DNA path - Removed the
flexiformattermodule. It existed to move the barcode and UMI from the read name into BAM tags, which is no longer needed: flexiplex writesCB/CR/UB/URinto the FASTQ comment itself, andminimap2 -ynow carries them onto every alignment.-yis applied only when demultiplexing with flexiplex, since a basecaller comment on any other FASTQ would be copied in verbatim and produce unparseable records - Deduplication now groups on
XB/UBrather thanCB/UR.XBis a derived tag holding the corrected barcode for a called cell and the uncorrected one for a droplet flexiplex could not match to the known list, so uncalled droplets stay separate instead of pooling into a single-cell.UBis flexiplex's corrected UMI;URholds the raw one. Theblazepath writes both tags too, so the two demultiplexers stay comparable. DNApicard MarkDuplicateskeys onXBfor the same reason - Reads in which flexiplex found no barcode region at all are now discarded during assignment rather than carried through the pipeline. They belong to no droplet, called or empty, and flexiplex gives them a one-character UMI placeholder instead of a real UMI. Dropping them at the FASTQ stage keeps them out of alignment entirely and means every alignment downstream has a real
CR, a realXBand a full-width UMI - Added a second IsoQuant run over the same alignments, published to
isoquant_all_droplets/, grouping onXBso that droplets below the knee are quantified alongside the called cells. The existingisoquant/run is unchanged and still groups onCB - Added the lr-kallisto equivalent of that second pass, published to
lrkallisto_all_droplets/: a newlrkallisto_countsubworkflow runs the counting half (split_bc_umi->kallisto bus->bustools->kallisto quant-tcc) twice off the same reads and the same index, once keyed on the barcode flexiplex matched to the known list and once onXB. The existinglrkallisto/outputs are unchanged - Bounded that all-droplet pass.
XBis the raw uncorrected barcode below the knee and the pass deliberately skipsbustools correct, so nothing limited the column count and sequencing error alone minted a barcode per read — on a full sample ~14.1M columns against ~8.8k called cells, which the long-read EM could not finish.split_bc_uminow takes an allowlist built from the read counts flexiplex has already collected, so barcodes below--lrkallisto_all_min_readsare dropped before pseudoalignment. At the default of 100 that is ~233k columns retaining ~54% of the below-knee reads, and it boundskallisto busandbustoolsas well as the EM - The lr-kallisto
gene/matrix is the integer UMI count matrix frombustools count --genecounts --umi-gene, in both thelrkallisto/andlrkallisto_all_droplets/passes. Thekallisto quant-tccEM gene abundances were considered for it and rejected: they are fractional (0.25, 0.333333, ...), which CellBender and every other tool with a count likelihood reject, and they inflatenFeaturein the Seurat QC because any nonzero EM weight counts as a detected gene. The all-droplet matrix exists for exactly that ambient/empty-droplet use. The published matrix is the standard kb-python gene matrix: a UMI compatible with more than one gene is discarded rather than distributed, so gene totals are lower than the EM totals by the gene-ambiguous fraction (28% on the chr21 test data, a worst case for so small a reference).--multimappingis deliberately not used for the gene count, since with--genecountsit splits ambiguous UMIs evenly and reintroduces the fractions; it can only be added through the separateext.args2ofBUSTOOLS_TCC, and the module then fails rather than publish a fractional matrix under anintegerheader.gene/features.tsv.gzcarries gene id, gene name and feature type, so the directory reads as a Cell Ranger v3 MEX directory. The EM gene abundances are still published asquant/matrix.abundance.gene.mtx. Thetranscript/matrix (and the oarfish matrix) remain EM estimates and stay fractional, so their Seurat QCnFeaturestill counts every nonzero EM weight and is not comparable with the IsoQuant transcript QC - The cDNA genome alignment is now produced whenever there are cDNA reads, rather than only when a genome quantifier was requested. Previously
--quantifier lrkallistoskipped the entire cDNA alignment arm, so an lr-kallisto run produced no cDNA BAM at all; the BAM is a deliverable in its own right (coverage, IGV, fusions, variants) and is independent of how the reads are quantified. lr-kallisto still pseudoaligns from the FASTQ and consumes nothing from this alignment
Parameter changes
--whitelistreplaced by--whitelist_dnaand--whitelist_cdna- Bundled whitelists in
assets/changed from.zipto.gz - Added
--split_amountparameter to control FASTQ/barcode splitting for parallelization - Added
--kallisto_kmer_size(default:63),--kallisto_threshold(default:0.8),--kallisto_platform(default:ONT) and--kallisto_dlist(default:true) for thelrkallistoquantifier - Added
--lrkallisto_all_min_reads(default:100), the minimum read count for a barcode to enter the lr-kallisto all-droplet matrix.0restores the previous unbounded behaviour - Added
--skip_lrkallisto_all_droplets(default:false) to skip the lr-kallisto all-droplet counting pass entirely - Added
--skip_cdna_dedup(default:false), which turns offumi_toolsdeduplication on the cDNA arm only.--skip_dedupcontinues to govern the DNA arm'spicard MarkDuplicates, so the now-always-on cDNA deduplication chain can be turned down without silently un-deduplicating the DNA BAM - Added
--dedup_per_gene(default:true, genome alignments) and--dedup_per_contig(default:true, transcriptome alignments) to control UMI grouping.--dedup_per_generequires--gtf
Fixes
- Fixed the three
flexiplexmodules declaringconda "${moduleDir}/environment.yml"when the file on disk isenvironment.yaml, which left them unresolvable under-profile conda split_gene_statusno longer emits an empty BAM.umi_toolsaborts on one, and a contig with no annotation produces exactly that — 50 of 85 contigs on a real GRCh38 run- Fixed
split_amountparameter type coercion so it is read as an integer under the Nextflow v2 strict syntax parser - Fixed the transcriptome deduplication path dropping the sample
typefrom the metadata when splitting the BAM by transcript group, which published the deduplicated BAM, its samtools stats and the oarfish outputs under<sample>/null/transcriptome/instead of<sample>/<type>/transcriptome/ - Corrected
docs/output.md, which listed the deduplicated BAMs underbam/dedup_umitools/andbam/dedup_picard/while the pipeline publishes them underbam/dedup/ - Fixed the
10X_5v3UMI width used by thelrkallistoquantifier (10, not 12), which otherwise caused every read to be dropped bysplit_bc_umi.py's length check - Fixed
SAMTOOLS_INDEX_DEDUPin thededup_umissubworkflow indexing theumi_toolsoutput channel unconditionally, which left the--dedup_tool picardbranch indexing a channel that was never populated - The
picard MarkDuplicatesduplication metrics are no longer discarded on the DNA branch.PICARD_MARKDUPLICATES.out.metricswas emitted by the module but published nowhere and never passed to MultiQC, so a DNA run reported no duplicate rate at all. The metrics are now published to<sample>/dna/qc/dedup/*.metrics.txtand summarised in the MultiQC report - Added
analysis/dna_dedup/, a standalone QC script and write-up quantifying what MarkDuplicates keys on for ONT single-end reads. It confirms the soft-clip problem behind the per-gene UMI grouping change does not affect the DNA branch, because the variable adapter clip sits at the read's 3' end while MarkDuplicates anchors on the 5' end --quantifiernow defaults tonullrather than an empty string.--quantifierhas been optional since DNA support was added, since DNA samples are aligned and deduplicated but never quantified, but the empty default was still checked against the parameter's own regex and failed it. Every DNA-only run therefore aborted at parameter validation with"" does not match regular expression,-profile test_dnaincluded. The samplesheet-aware check that errors only when acdnarow is present and no quantifier was given is now what actually reports the problem. Addedtests/dna.nf.testso the DNA-only path is covered by nf-test- Removed four module includes from the
process_longread_scrnasubworkflow that were never invoked (PICARD_MARKDUPLICATES,SAMTOOLS_FLAGSTAT_DEDUP,SAMTOOLS_INDEX_DEDUP,SAMTOOLS_FILTER_DEDUP), along with theconf/modules.configselector forSAMTOOLS_FILTER_DEDUP, which matched no process. The underlying modules are all still used elsewhere; only the dead aliases are gone - DNA samtools QC files are no longer named with a doubled type suffix. The
BAM_STATS_SAMTOOLSprefix appended a literal.dnaon top of themeta.typeinterpolation, so a DNA sample produced<id>.dna.dna.flagstat. They are now<id>.dna.flagstat/.idxstats/.stats. File contents are unchanged; MultiQC sample labels lose the duplicate suffix too. ThePICARD_MARKDUPLICATESprefix carried the same doubled literal and is fixed alongside it, so the newly published duplication metrics are<id>.dna.dedup.metrics.txt - Fixed two
tests/.nftignorepatterns that had silently stopped matching.**/genome/isoquant/*.mtxand itsisoquant_all_dropletscounterpart were added to keep the count matrices out of content checking, but the published path later gained two nested<sample>.<type>/levels, so the matrices were being md5-compared again. Also added the DNA dedup BAM and its index, which are not byte-reproducible between runs --samtools sortties are resolved differently across runs -- and whose cDNA equivalents were already ignored. Both files are still checked by name. The newly published MarkDuplicates metrics are ignored for content on the same grounds: the file carries a# Started on:timestamp and the full work-directory command line, which is why the nf-core module test snapshots only its first three lines - Fixed flexiplex being unable to demultiplex
10X_5v3and10X_3v4cDNA. The explicit 5' v3 pattern was guarded bybarcode_format == "10x5v3", which never matches the schema value10X_5v3, so both formats fell through to a derived-d 10x5v3/-d 10x3v4preset that flexiplex does not have. Every format is now mapped explicitly:10X_3v4uses-d 10x3v3(same geometry, different whitelist) and10X_5v3uses the 5' flanks with a 10 bp UMI - Transcriptome
umi_tools dedupnow gets--per-gene --per-contiginstead of--per-contigalone, which umi_tools rejects (need to use --per-gene with --per-contig). Every run with a transcriptome quantifier (oarfish) failed at this step. No test profile exercises oarfish, so it was found by the real-data release test
Known issues
- The called-cell IsoQuant (
isoquant/) and oarfish matrices contain a-barcode column. Running flexiplex with-a trueleavesCBas-on reads that matched no known barcode, and these passes group onCB, so all such reads pool into one pseudo-cell. On the release test (one 10X Multiome cDNA library, 4,379 called cells) it held 51-55% of all counts. Drop the-column before downstream analysis. Thelrkallisto/matrices and the all-droplet passes (isoquant_all_droplets/,lrkallisto_all_droplets/, keyed onXB) are not affected group_transcripts.shtakes the transcript id from GTF field 12. In a versioned Ensembl/Cell Ranger GTF that field is thegene_versionvalue, so every transcript is assigned tochr_Unknownand transcriptome deduplication runs as a single task instead of one per chromosome. Results are unaffected; only parallelism is lost. GENCODE GTFs, wheretranscript_idis the second attribute, are grouped correctly. This behaviour is inherited from nf-core/scnanoseq- No nf-test profile runs
oarfish, so the transcriptome deduplication and oarfish path is covered only by-profile test_full OARFISHcarries theprocess_lowlabel (12 GB per attempt) but peaked at 36.3 GB on the release test. It finished only because the institutional profile allowed further retries; with the pipeline's ownmaxRetries = 1it stops at 24 GB and fails. Until the label is raised, give it more memory with-c, e.g.process { withName: '.*:OARFISH' { memory = 48.GB } }