Skip to content

Releases: IntGenomicsLab/lrsinglecell

Release list

v1.0.0

Choose a tag to compare

@ljwharbers ljwharbers released this 31 Oct 15:16
0a643f5

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 is INTGENOMICSLAB_LRSINGLECELL and the software versions file is pipeline_info/lrsinglecell_software_mqc_versions.yml
  • Test profiles read their data from IntGenomicsLab/test-datasets (lrsinglecell/), copied unchanged from nf-core/test-datasets scnanoseq
  • Moved to the nf-core/tools 4.2.1 template (minimum Nextflow 25.10.4, nf-schema 2.7.3, nft-utils 1.2.1, updated utils_* subworkflows, multiqc module and fastqc module (FastQC 0.13.0, only used with --skip_fastqc false))

Enhancements

  • DNA samples are now deduplicated by mark_dna_duplicates.py instead of picard 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; on chrM/MT only 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; see docs/output.md
  • Added support for processing single-cell/nuclei DNA samples alongside cDNA samples in the same pipeline run, via a new optional type column (dna/cdna) in the input samplesheet
  • Added flexiplex/discovery, flexiplex/filter and flexiplex/assign modules for barcode extraction, filtering and assignment
  • Added demultiplex_flexiplex and demultiplex_blaze subworkflows, replacing the previous inline demultiplexing steps in the main workflow
  • Added align_deduplicate_dna subworkflow (minimap2 alignment, mark_dna_duplicates.py, BAM_STATS_SAMTOOLS) for DNA sample processing
  • Added demux_tool_cdna/demux_tool_dna parameters to select flexiplex or blaze for cDNA demultiplexing (DNA currently supports flexiplex only)
  • Added lrkallisto as a third --quantifier option: an alignment-free path that pseudoaligns the demultiplexed FASTQ with lr-kallisto (kallisto bus --long), counts genes with bustools count --genecounts and quantifies transcripts with kallisto quant-tcc --long, deduplicating UMIs with bustools in both cases. Adds a test_lrkallisto profile and the gffread module 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_genes and split_gene_status modules, GX/GN/GS tags) and grouped with umi_tools --per-gene; reads with an ambiguous gene call or none are deduplicated positionally and merged back, so no reads are dropped. Transcriptome alignments use umi_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-chunk umi_tools dedup runs by the new umitools_dedup_summary module, and shown as a "UMI-tools Deduplication" table in the MultiQC report. Previously the umi_tools logs were discarded and no duplication rate was reported
  • umi_tools dedup now always takes the cell barcode from the corrected CB tag. Previously this applied only when --demux_tool_cdna flexiplex was set; the blaze path 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 in CR, 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 flexiformatter module. It existed to move the barcode and UMI from the read name into BAM tags, which is no longer needed: flexiplex writes CB/CR/UB/UR into the FASTQ comment itself, and minimap2 -y now carries them onto every alignment. -y is 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/UB rather than CB/UR. XB is 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. UB is flexiplex's corrected UMI; UR holds the raw one. The blaze path writes both tags too, so the two demultiplexers stay comparable. DNA picard MarkDuplicates keys on XB for 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 real XB and a full-width UMI
  • Added a second IsoQuant run over the same alignments, published to isoquant_all_droplets/, grouping on XB so that droplets below the knee are quantified alongside the called cells. The existing isoquant/ run is unchanged and still groups on CB
  • Added the lr-kallisto equivalent of that second pass, published to lrkallisto_all_droplets/: a new lrkallisto_count subworkflow 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 on XB. The existing lrkallisto/ outputs are unchanged
  • Bounded that all-droplet pass. XB is the raw uncorrected barcode below the knee and the pass deliberately skips bustools 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_umi now takes an allowlist built from the read counts flexiplex has already collected, so barcodes below --lrkallisto_all_min_reads are dropped before pseudoalignment. At the default of 100 that is ~233k columns retaining ~54% of the below-knee reads, and it bounds kallisto bus and bustools as well as the EM
  • The lr-kallisto gene/ matrix is the integer UMI count matrix from bustools count --genecounts --umi-gene, in both the lrkallisto/ and lrkallisto_all_droplets/ passes. The kallisto quant-tcc EM 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 inflate nFeature in 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). --multimapping is deliberately not used for the gene count, since with --genecounts it splits ambiguous UMIs evenly and reintroduces the fractions; it can only be added through the separate ext.args2 of BUSTOOLS_TCC, and the module then fails rather than publish a fractional matrix under an integer header. gene/features.tsv.gz carries 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 as quant/matrix.abundance.gene.mtx. The transcript/ matrix (and the oarfish matrix) remain EM estimates and stay fractional, so their Seurat QC nFeature still 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 lrkallisto skipped 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 th...
Read more