diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index bccfa04..cf1b432 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -36,7 +36,6 @@ jobs: shell: bash -l {0} run: | pip install . - pip install "pulp==2.7.0" --no-deps # pin pulp to avoid dependency conflicts - name: Install dependencies shell: bash -l {0} diff --git a/README.rst b/README.rst index c612c90..6706822 100644 --- a/README.rst +++ b/README.rst @@ -8,9 +8,9 @@ .. image:: https://coveralls.io/repos/github/sequana/demultiplex/badge.svg?branch=main :target: https://coveralls.io/github/sequana/demultiplex?branch=main -.. image:: https://img.shields.io/badge/python-3.8%20%7C%203.9%20%7C3.10-blue.svg +.. image:: https://img.shields.io/badge/python-3.9%20%7C3.10%20%7C3.11-blue.svg :target: https://pypi.python.org/pypi/sequana - :alt: Python 3.8 | 3.9 | 3.10 + :alt: Python 3.9 | 3.10 | 3.11 .. image:: http://joss.theoj.org/papers/10.21105/joss.00352/status.svg :target: http://joss.theoj.org/papers/10.21105/joss.00352 @@ -116,6 +116,7 @@ Changelog ========= ======================================================================= Version Description ========= ======================================================================= +1.6.0 * implement bclconvert along bcl2fastq (use --use-bcl-convert) 1.5.2 * rename requirements.txt into tools.txt; update __init__ 1.5.1 * switch working directory to fastq instead of demultiplex(regression) 1.5.0 * Uses click and new sequana_pipetools, add multiqc diff --git a/pyproject.toml b/pyproject.toml index 06dc7f1..2176527 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "poetry.core.masonry.api" [tool.poetry] name = "sequana-demultiplex" -version = "1.5.3" +version = "1.6.0" description = "Pipeline that runs bcl2fastq and ease demultiplexing of Sequencing data" authors = ["Sequana Team"] license = "BSD-3" @@ -19,7 +19,6 @@ classifiers = [ "Intended Audience :: Science/Research", "License :: OSI Approved :: BSD License", "Operating System :: POSIX :: Linux", - "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", diff --git a/sequana_pipelines/demultiplex/config.yaml b/sequana_pipelines/demultiplex/config.yaml index 92fe814..da006d4 100644 --- a/sequana_pipelines/demultiplex/config.yaml +++ b/sequana_pipelines/demultiplex/config.yaml @@ -31,6 +31,8 @@ general: apptainers: graphviz: "https://zenodo.org/record/7928262/files/graphviz_7.0.5.img" + bcl2fastq: "https://zenodo.org/record/13304101/files/bcl2fastq_2.20.0.img" + bcl_convert: "https://zenodo.org/record/10934031/files/bclconvert_4.2.7.img" cellranger_atac: "https://zenodo.org/record/8423332/files/cellranger_atac_2.1.0.img" multiqc: https://zenodo.org/record/10205070/files/multiqc_1.16.0.img diff --git a/sequana_pipelines/demultiplex/demultiplex.rules b/sequana_pipelines/demultiplex/demultiplex.rules index a9a0ed5..979ce1f 100644 --- a/sequana_pipelines/demultiplex/demultiplex.rules +++ b/sequana_pipelines/demultiplex/demultiplex.rules @@ -52,6 +52,8 @@ if config["general"]["mode"] == "bcl2fastq": shell(f"sequana samplesheet --check {samplesheet}") +mode = config["general"]["mode"] + rule pipeline: input: "Stats/summary.txt", @@ -88,7 +90,7 @@ rule plot_unknown_barcodes: -if config["general"]["mode"] == "bcl2fastq": +if mode == "bcl2fastq": rule bcl2fastq: input: @@ -108,10 +110,12 @@ if config["general"]["mode"] == "bcl2fastq": threads: config['bcl2fastq']['threads'] resources: **config['bcl2fastq']['resources'] + container: + config['apptainers']['bcl2fastq'] shell: manager.get_shell("bcl2fastq/run", "v1") -elif config["general"]["mode"] == "cellranger_atac": +elif mode == "cellranger_atac": rule cellranger_atac: input: @@ -123,7 +127,7 @@ elif config["general"]["mode"] == "cellranger_atac": outdir=".", indir= config["general"]["input_directory"], container: - "https://zenodo.org/record/8423332/files/cellranger_atac_2.1.0.img" + config['apptainers']['cellranger_atac'] resources: **config['cellranger_atac']['resources'] shell: @@ -134,45 +138,86 @@ elif config["general"]["mode"] == "cellranger_atac": """ -elif config["general"]["mode"] == "bclconvert": +elif mode == "bclconvert": rule bclconvert: input: samplesheet=config['general']['samplesheet_file'] output: - "Stats/Stats.json", + "Reports/Demultiplex_Stats.csv", params: - indir= config["general"]["input_directory"], - #ignore_missing_bcls=config['bclconvert']['ignore_missing_bcls'], - bcl_sampleproject_subdirectories=config['bclconvert']['bcl_sampleproject_subdirectories'], - merge_all_lanes=config['bclconvert']['merge_all_lanes'], + indir=config["general"]["input_directory"], + bcl_sampleproject_subdirectories="true" if config['bclconvert']['bcl_sampleproject_subdirectories'] else "false", + no_lane_splitting="--no-lane-splitting true" if config['bclconvert']['merge_all_lanes'] else "--no-lane-splitting false", options=config['bclconvert']['options'] threads: config['bclconvert']['threads'] resources: **config['bclconvert']['resources'] + container: + config['apptainers']['bcl_convert'] shell: """ - - PWD = `pwd` - bcl-convert --sample-sheet {input.samplesheet} --bcl-input-directory {params.indir} \ - --output-directory ${PWD} --force \ - # make sure to use those options for cluster usage. indeed default is number of available CPU - # - --bcl-num-conversion-threads {snakemake.threads} \ - --bcl-num-compression-threads {snakemake.threads} \ - # command line --ignore-missing-bcls not supported anymore. Must be in the sample sheet as - # BarcodeMismatchIndex1,# and BarcodeMismatchIndex2,# - --no-lane-splitting {params.merge_all_lanes} - # " --no-bgzf-compression " does not exist anymore - # - # unlike bcl2fastq, sample project is ignore and all fastq stored in 1 main directory. - # use --bcl-sampleproject-subdirectories we have same behaviour as in bcl2fastq - --bcl-sampleproject-subdirectories {params.bcl_sampleproject_subdirectories} - # SAMPLE ID is used to name the files ; sample name is not used except if used with - # --sample-name-column-enabled true in which case the fastq files are - # name project/ID/Name.fastq.gz - + bcl-convert --sample-sheet {input.samplesheet} \ + --bcl-input-directory {params.indir} \ + --output-directory . \ + --force \ + --bcl-num-conversion-threads {threads} \ + --bcl-num-compression-threads {threads} \ + --bcl-sampleproject-subdirectories {params.bcl_sampleproject_subdirectories} \ + {params.no_lane_splitting} \ + {params.options} """ + rule bclconvert_to_stats_json: + """Convert bclconvert Reports CSVs into the Stats/Stats.json format + expected by the sequana StatsFile class, so all downstream plot rules + can be reused without modification.""" + input: + demux="Reports/Demultiplex_Stats.csv", + output: + "Stats/Stats.json", + run: + import json, os + import pandas as pd + + demux = pd.read_csv(input.demux) + + barcodes_file = "Reports/Top_Unknown_Barcodes.csv" + barcodes = pd.read_csv(barcodes_file) if os.path.exists(barcodes_file) else pd.DataFrame(columns=["Lane", "index", "index2", "# Reads"]) + + # Build ConversionResults grouped by lane + conversion_results = [] + for lane_num, lane_df in demux.groupby("Lane"): + demux_results = [] + undetermined_reads = 0 + for _, row in lane_df.iterrows(): + if row["SampleID"] == "Undetermined": + undetermined_reads = int(row["# Reads"]) + else: + demux_results.append({ + "SampleId": row["SampleID"], + "NumberReads": int(row["# Reads"]), + }) + conversion_results.append({ + "LaneNumber": int(lane_num), + "DemuxResults": demux_results, + "Undetermined": {"NumberReads": undetermined_reads}, + }) + + # Build UnknownBarcodes grouped by lane + unknown_barcodes = [] + for lane_num, lane_df in barcodes.groupby("Lane"): + barcode_dict = {} + for _, row in lane_df.iterrows(): + idx2 = row.get("index2", "") + key = f"{row['index']}+{idx2}" if pd.notna(idx2) and str(idx2).strip() else str(row["index"]) + barcode_dict[key] = int(row["# Reads"]) + unknown_barcodes.append({"Lane": int(lane_num), "Barcodes": barcode_dict}) + + os.makedirs("Stats", exist_ok=True) + with open(output[0], "w") as fh: + json.dump({"ConversionResults": conversion_results, + "UnknownBarcodes": unknown_barcodes}, fh, indent=2) + rule plot_barplot_samples: input: "Stats/Stats.json" @@ -197,23 +242,28 @@ rule plot_summary: s.to_summary_reads(output.summary) - # ========================================================== multiqc multiqc_params_options = config['multiqc']['options'] if manager.config.multiqc.config_file: multiqc_params_options += f" -c {manager.config.multiqc.config_file}" +if mode == "bclconvert": + _multiqc_input = "Reports/Demultiplex_Stats.csv" + _multiqc_module = "bclconvert" +else: + _multiqc_input = ["Stats/Stats.json", "Stats/summary.txt"] + _multiqc_module = config['multiqc']['modules'] + rule multiqc: input: - "Stats/Stats.json", - "Stats/summary.txt", + _multiqc_input, output: "multiqc/multiqc_report.html" params: options=multiqc_params_options, input_directory=config['multiqc']['input_directory'], config_file=config['multiqc']['config_file'], - modules=config['multiqc']['modules'] + modules=_multiqc_module log: "multiqc/multiqc.log" resources: @@ -262,7 +312,7 @@ onsuccess: image1 = SequanaReport.png_to_embedded_png("dummy", "barcodes.png", - style="text-align:center; width:60%; height:40%", + style="text-align:center; width:60%; height:40%", alt="barcodes") image2 = SequanaReport.png_to_embedded_png("dummy", "summary.png", @@ -271,24 +321,35 @@ onsuccess: image3 = SequanaReport.png_to_embedded_png("dummy", "samples.png", style="width:60%; height:40%", alt="sample") - intro = """
-

This report summarizes the demultiplexing of your raw data. Some additional help to interpret the following plots can be found on the pipeline home page and wiki.

+

Overview

+

This report summarizes the demultiplexing of your raw BCL data into per-sample FASTQ files. +The FASTQ files are located in the current working directory, organised by project and sample as defined in your sample sheet. +Undetermined reads (those that could not be assigned to any sample) are written to files prefixed with Undetermined. +A detailed MultiQC report is available here. +Additional help on interpreting the results can be found on the pipeline +home page and wiki.

- -

The following image shows the indices found in the most undetermined barcodes per lane. An index in excess may indicate a wrong sample sheet with a potential typo or simply a missing sample/index. -The 20 most represented undetermined barcodes are available here. A full list of all undetermined indices is also available here as a CSV file. The vertical dashed red line indicates the median number of reads found in the determined samples. If the bars are well below this line, your data is probably well demultiplexed.

+

Top undetermined barcodes per lane

+

The plot below shows the most frequent index sequences found among the undetermined reads, broken down by lane. +A barcode appearing in large numbers may indicate a typo in the sample sheet, a missing sample, or an index that was not declared. +The vertical dashed red line marks the median read count across all determined samples: bars significantly below this line suggest the run is well demultiplexed, while bars approaching or exceeding it warrant investigation. +The top 20 undetermined barcodes are available as a CSV file here, and the complete list of all undetermined indices is available here.


{}
- -

The following image shows the ratio of determined vs. undetermined reads after demultiplexing. An excess of undetermined reads (more than 10-20%) may indicate a wrongly labelled sample in your sample sheet.

+

Determined vs. undetermined reads

+

The plot below shows the proportion of reads that were successfully assigned to a sample (determined) versus those that could not be assigned (undetermined), for each lane. +A high fraction of undetermined reads (typically above 10–20%) may indicate a mislabelled sample, an index swap, or an error in the sample sheet. +In a well-demultiplexed run, the vast majority of reads should fall in the determined category.


{}
- -

The following image shows the number of reads per sample.

+

Read counts per sample

+

The plot below shows the total number of reads assigned to each sample after demultiplexing. +Samples with unexpectedly low read counts may have been affected by index mismatches or low cluster density on the corresponding lane. +Large imbalances between samples loaded at equal proportions may indicate pipetting errors during library preparation.


{}
""".format(image1, image2, image3) diff --git a/sequana_pipelines/demultiplex/main.py b/sequana_pipelines/demultiplex/main.py index b0cf944..77929c6 100644 --- a/sequana_pipelines/demultiplex/main.py +++ b/sequana_pipelines/demultiplex/main.py @@ -27,7 +27,19 @@ help = init_click( NAME, groups={ - "Pipeline Specific": ["--method", "--skip-multiqc"], + "Pipeline Specific": [ + "--threads", + "--barcode-mismatch", + "--merging-strategy", + "--bcl-directory", + "--sample-sheet", + "--use-bcl-convert", + "--no-ignore-missing-bcls", + "--strict-mode", + "--bgzf-compression", + "--mars-seq", + "--scatac-seq", + ], }, ) @@ -35,7 +47,6 @@ @click.command(context_settings=help) @include_options_from(ClickSnakemakeOptions, working_directory="fastq") @include_options_from(ClickSlurmOptions) -@include_options_from(ClickInputOptions, add_input_readtag=False) @include_options_from(ClickGeneralOptions) @click.option( "--threads", @@ -102,6 +113,14 @@ bcl2fastq uses this option; By default we don't. Set --bgzl--compression flag to set it back. For bcl2fastq only.""", ) +@click.option( + "--use-bcl-convert", + "use_bcl_convert", + is_flag=True, + default=False, + show_default=True, + help="Use bcl-convert instead of bcl2fastq for demultiplexing.", +) @click.option( "--mars-seq", default=False, @@ -181,12 +200,10 @@ def main(**options): cfg.bclconvert.threads = options.threads cfg.bcl2fastq.strict_mode = not options.strict_mode - if options.merging_strategy == "merge": - cfg.bcl2fastq.merge_all_lanes = True - elif options.merging_strategy in ["none", "none_and_force"]: - cfg.bcl2fastq.merge_all_lanes = False + merge = options.merging_strategy == "merge" + cfg.bcl2fastq.merge_all_lanes = merge + cfg.bclconvert.merge_all_lanes = merge - # if options.mars_seq: cfg.bcl2fastq.options = " --minimum-trimmed-read-length 15 --mask-short-adapter-reads 15 " if options.merging_strategy in ["merge"]: @@ -196,9 +213,17 @@ def main(**options): elif options.scatac_seq: cfg.cellranger_atac.options = "" cfg.general.mode = "cellranger_atac" - #elif options.bclconvert: # All other cases with bcl2fastq - # cfg.general.mode = "bclconvert" - + elif options.use_bcl_convert: + cfg.bclconvert.threads = options.threads + cfg.bclconvert.strict_mode = options.strict_mode + cfg.general.mode = "bclconvert" + # bcl-convert writes logs to /var/log/bcl-convert which is read-only inside + # Apptainer/Singularity containers. Bind a local writable directory instead. + if options.apptainer_prefix: + log_dir = os.path.abspath("bcl_convert_logs") + os.makedirs(log_dir, exist_ok=True) + bind = f"--bind {log_dir}:/var/log/bcl-convert" + options.apptainer_args = f"{options.apptainer_args} {bind}".strip() else: cfg.general.mode = "bcl2fastq" try: diff --git a/test/data/SampleSheet_bclconvert.csv b/test/data/SampleSheet_bclconvert.csv new file mode 100644 index 0000000..e9c3590 --- /dev/null +++ b/test/data/SampleSheet_bclconvert.csv @@ -0,0 +1,18 @@ +[Header] +FileFormatVersion,2 +RunName,Test +InstrumentPlatform,NextSeq/MiniSeq +Description,TEST + +[Reads] +Read1Cycles,36 + +[BCLConvert_Settings] +SoftwareVersion,4.2.7 +AdapterRead1,AGATCGGAAGAGCACACGTCTGAACTCCAGTCA +AdapterRead2,AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT + +[BCLConvert_Data] +Sample_ID,index,Sample_Project +Sample1,ATCACG,TEST +Sample2,CGATGT,TEST diff --git a/test/test_main.py b/test/test_main.py index 57cdceb..3e6a664 100644 --- a/test/test_main.py +++ b/test/test_main.py @@ -10,6 +10,7 @@ from . import test_dir samplesheet = f"{test_dir}/data/SampleSheet.csv" +samplesheet_bclconvert = f"{test_dir}/data/SampleSheet_bclconvert.csv" bcldir = f"{test_dir}/data/bcl" def test_help(): @@ -37,6 +38,17 @@ def test_standalone_script(tmp_path): assert results.exit_code == 0 +def test_standalone_bclconvert(tmp_path): + directory = tmp_path / "test" + directory.mkdir() + + runner = CliRunner() + results = runner.invoke(main, ["--bcl-directory", bcldir, "--merging-strategy", "merge", + "--working-directory", str(directory), "--force", + "--sample-sheet", samplesheet_bclconvert, "--use-bcl-convert"]) + assert results.exit_code == 0 + + def test_standalone_baddies(tmp_path): directory = tmp_path / "test" directory.mkdir()