From bb9acf55bea4d1ce29bb4d85960dfe35d41cb801 Mon Sep 17 00:00:00 2001 From: Neil Horner Date: Tue, 13 Feb 2024 08:48:21 +0000 Subject: [PATCH] Resolve CW-3468 "Incorrect ids" --- .pre-commit-config.yaml | 2 +- CHANGELOG.md | 2 +- README.md | 7 ++ bin/workflow_glue/report.py | 109 ++++++---------------------- docs/06_input_parameters.md | 7 ++ main.nf | 71 ++++++++---------- nextflow.config | 2 +- subworkflows/JAFFAL/gene_fusions.nf | 2 +- subworkflows/reference_assembly.nf | 5 ++ 9 files changed, 79 insertions(+), 128 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 2ca2e89..2f0bb9f 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -8,7 +8,7 @@ repos: always_run: true pass_filenames: false additional_dependencies: - - epi2melabs>=0.0.51 + - epi2melabs>=0.0.52 - id: build_models name: build_models entry: datamodel-codegen --strict-nullable --base-class workflow_glue.results_schema_helpers.BaseModel --use-schema-description --disable-timestamp --input results_schema.yml --input-file-type openapi --output bin/workflow_glue/results_schema.py diff --git a/CHANGELOG.md b/CHANGELOG.md index d79b409..d77e00f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,7 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [unreleased] +## [v1.1.0] ### Changed - Improve documentation around filtering of transcripts done before DTU analysis. - Renamed files: diff --git a/README.md b/README.md index 8b1b42e..b96f834 100644 --- a/README.md +++ b/README.md @@ -122,6 +122,13 @@ input_reads.fastq ─── input_directory ─── input_directory | analyse_unclassified | boolean | Analyse unclassified reads from input directory. By default the workflow will not process reads in the unclassified directory. | If selected and if the input is a multiplex directory the workflow will also process the unclassified directory. | False | +### Output Options + +| Nextflow parameter name | Type | Description | Help | Default | +|--------------------------|------|-------------|------|---------| +| out_dir | string | Directory for output of all user-facing files. | | output | + + ### Sample Options | Nextflow parameter name | Type | Description | Help | Default | diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py index cdc7afc..787161c 100755 --- a/bin/workflow_glue/report.py +++ b/bin/workflow_glue/report.py @@ -51,23 +51,17 @@ def argparser(): "--gff_annotation", required=False, nargs='+', help="transcriptome annotation gff file") parser.add_argument( - "--gffcompare_dir", required=False, default=None, nargs='*', + "--gffcompare_dir", required=False, default=None, help="gffcompare outout dir") parser.add_argument( "--pychop_report", required=False, default=None, help="TSV summary file of pychopper statistics") - parser.add_argument( - "--sample_ids", required=True, nargs='+', - help="List of sample ids") parser.add_argument( "--isoform_table", required=False, type=Path, help="Path to directory of TSV files with isoform summaries") parser.add_argument( "--isoform_table_nrows", required=False, type=int, default=5000, help="Maximum rows to display in isoforms table") - parser.add_argument( - "--cluster_qc_dirs", required=False, type=str, default=None, nargs='*', - help="Directory with various cluster quality csvs") parser.add_argument( "--jaffal_csv", required=False, type=str, default=None, help="Path to JAFFAL results csv") @@ -345,7 +339,7 @@ def grouped_bar(df, title="", tilted_xlabs=False): return p -def gff_compare_plots(report, gffcompare_outdirs, sample_ids): +def gff_compare_plots(report, gffcompare_outdirs): """Create various sections and plots in a WfReport. :param report: aplanat WFReport @@ -383,7 +377,10 @@ def gff_compare_plots(report, gffcompare_outdirs, sample_ids): tabs = [] gff_fails = False - for id_, dir_ in zip(sample_ids, gffcompare_outdirs): + sample_ids = [] + for dir_ in gffcompare_outdirs: + sample_id = dir_.name + sample_ids.append(sample_id) # Get sample ids fromt the folder name stats, _, miss, novel, total = \ parse_gffcmp_stats(dir_ / 'str_merged.stats') @@ -396,7 +393,7 @@ def gff_compare_plots(report, gffcompare_outdirs, sample_ids): tabs.append(Panel( child=gridplot( [bar_totals, bar_performance, bar_missed, bar_novel], - ncols=2, width=350, height=260), title=id_)) + ncols=2, width=350, height=260), title=sample_id)) else: gff_fails = True @@ -445,7 +442,7 @@ def gff_compare_plots(report, gffcompare_outdirs, sample_ids): track_files = [x / 'str_merged.tracking' for x in gffcompare_outdirs] - df_tracking = load_sample_data( + df_tracking = load_data_add_sample_id( track_files, sample_ids, read_func=lambda x: pd.read_csv( x, sep="\t", header=None, @@ -531,7 +528,7 @@ def gff_compare_plots(report, gffcompare_outdirs, sample_ids): sys.stderr("Cannot find .tmap files in {}".format(gffcompare_outdirs)) return - df_tmap = load_sample_data(tmap_files, sample_ids) + df_tmap = load_data_add_sample_id(tmap_files, sample_ids) for id_, df in df_tmap.groupby('sample_id'): @@ -596,59 +593,6 @@ def pychopper_plots(report, pychop_report): section.plot(grid) -def cluster_quality(cluster_qc_dir, report, sample_ids): - """Make cluster quality section.""" - section = report.add_section() - section.markdown(''' - ### De novo clustering quality - - This section shows plots relating to the clustering quality performed - by isONclust2. The full length reads are mapped to a reference genome - to create a ground truth of reads mapped to clusters. This is then compared - to the de novo-generated clusters, and the following statistics are - generated. - - * [Homogeneity](https://scikit-learn.org/stable/modules/generated/ - sklearn.metrics.homogeneity_score.html): Penalises over-clustering. - - * [Completeness](https://scikit-learn.org/stable/modules/generated/ - sklearn.metrics.completeness_score.html): Penalises under-clustering. - - * [V-measure](https://clusteringjl.readthedocs.io/en/latest/vmeasure.html): - The harmonic mean of the homogeneity and completeness - - * [Adjusted Rand Index](https://scikit-learn.org/stable/modules/generated/ - sklearn.metrics.adjusted_rand_score.html): Intuitively, measures the - percentage of read pairs correctly clustered, normalized so that a perfect - clustering = 1 and a random cluster assignment achieves = 0 - - * NonSingleton: Clusters with multiple reads - * Singleton: Clusters consisting of a single read (These do not contribute - to the final transcript calling - I need to check this!) - - ''') - - tabs = [] - for id_, cluster_dir in zip(sample_ids, cluster_qc_dir): - plots = [] - for fn in ['v_ari_com_hom.csv', 'sing_nonsing.csv']: - # Skip the next two plots for now - # 'class_sizes1.csv', 'class_sizes2.csv']: - df = pd.read_csv(Path(cluster_dir) / fn) - bar = bars.simple_bar( - df.Statistic.values.tolist(), df.Value.values.tolist(), - colors=Colors.cerulean - ) - bar.xaxis.major_label_orientation = math.pi / 2.8 - plots.append(bar) - tabs.append(Panel( - child=gridplot(plots, ncols=4, - width=300, height=300), title=id_)) - - cover_panel = Tabs(tabs=tabs) - section.plot(cover_panel) - - def transcript_table(report, isoform_table, max_rows): """Create searchable table of transcripts. @@ -685,7 +629,7 @@ def transcript_table(report, isoform_table, max_rows): section.table(df, index=False) -def transcriptome_summary(report, gffs, sample_ids): +def transcriptome_summary(report, gffs): """ Plot transcriptome summaries. @@ -694,7 +638,6 @@ def transcriptome_summary(report, gffs, sample_ids): :param report: aplanat WFReport :param gffs: list of paths to gff transcriptome annotations - :param sample_ids: list of sample ids """ # test.db gets written to the git repo. section = report.add_section() @@ -703,7 +646,8 @@ def transcriptome_summary(report, gffs, sample_ids): ''') tabs = [] - for id_, gff in zip(sample_ids, gffs): + for gff in gffs: + sample_id = Path(gff).name plots = [] @@ -747,7 +691,7 @@ def transcriptome_summary(report, gffs, sample_ids): plots.append(bar_isos) box = bars.boxplot_series( - [id_] * len(transcript_lens), transcript_lens, + [sample_id] * len(transcript_lens), transcript_lens, width=70, ylim=(min(transcript_lens), max(transcript_lens)), title='transcript lengths') plots.append(box) @@ -780,14 +724,14 @@ def transcriptome_summary(report, gffs, sample_ids): tabs.append(Panel( child=gridplot(plots, ncols=4, - width=300, height=300), title=id_)) + width=300, height=300), title=sample_id)) cover_panel = Tabs(tabs=tabs) section.plot(cover_panel) -def load_sample_data(files, sample_ids, read_func=None): - """Load CSVs into dataframe, and assign sample_id column.""" +def load_data_add_sample_id(files, sample_ids, read_func=None): + """Load CSVs and concat into single into dataframe, and assign sample_id column.""" df_ = pd.DataFrame() if not files: return None @@ -890,18 +834,18 @@ def de_section(report): def main(args): """Run the entry point.""" - sample_ids = args.sample_ids - sample_ids.sort() - report = WFReport( "Transcript isoform report", "wf-transcriptomes", revision=args.revision, commit=args.commit) - # QC seq_stats_tabs(report, args.stats) if args.alignment_stats is not None: - df_aln_stats = load_sample_data(args.alignment_stats, sample_ids) + stats_dfs = [] + for stats_file in args.alignment_stats: + df = pd.read_csv(stats_file, sep='\t+') + stats_dfs.append(df) + aln_stats_df = pd.concat(stats_dfs) section = report.add_section() section.markdown(''' ### Read mapping summary @@ -910,28 +854,23 @@ def main(args): [seqkit](https://bioinf.shenwei.me/seqkit/) `seqkit bam -s`''') - section.table(df_aln_stats) + section.table(aln_stats_df) if args.pychop_report is not None: pychopper_plots(report, args.pychop_report) # Results if args.gff_annotation is not None: - transcriptome_summary( - report, args.gff_annotation, sample_ids) + transcriptome_summary(report, args.gff_annotation) if args.gffcompare_dir is not None: gff_compare_plots( report, - [Path(x) for x in args.gffcompare_dir], - sample_ids) + [x for x in Path(args.gffcompare_dir).iterdir()]) if args.isoform_table is not None: transcript_table(report, args.isoform_table, args.isoform_table_nrows) - if args.cluster_qc_dirs is not None: - cluster_quality(args.cluster_qc_dirs, report, sample_ids) - if args.de_report: de_section(report) diff --git a/docs/06_input_parameters.md b/docs/06_input_parameters.md index 22b58a2..e3a0efb 100644 --- a/docs/06_input_parameters.md +++ b/docs/06_input_parameters.md @@ -11,6 +11,13 @@ | analyse_unclassified | boolean | Analyse unclassified reads from input directory. By default the workflow will not process reads in the unclassified directory. | If selected and if the input is a multiplex directory the workflow will also process the unclassified directory. | False | +### Output Options + +| Nextflow parameter name | Type | Description | Help | Default | +|--------------------------|------|-------------|------|---------| +| out_dir | string | Directory for output of all user-facing files. | | output | + + ### Sample Options | Nextflow parameter name | Type | Description | Help | Default | diff --git a/main.nf b/main.nf index 5282271..b26246d 100644 --- a/main.nf +++ b/main.nf @@ -150,19 +150,18 @@ process preprocess_reads { tuple val("${meta.alias}"), path("${meta.alias}_pychopper_output/"), emit: pychopper_output - path "${meta.alias}_pychopper_output/${meta.alias}_pychopper.tsv", + path("${meta.alias}_pychopper_output/pychopper.tsv"), emit: report script: def cdna_kit = params.cdna_kit.split("-")[-1] - def extra_params = params.pychopper_opts ?: '' + def extra_params = params.pychopper_opts ?: '' """ pychopper -t ${params.threads} -k ${cdna_kit} -m ${params.pychopper_backend} ${extra_params} 'seqs.fastq.gz' ${meta.alias}_full_length_reads.fastq - mv pychopper.tsv ${meta.alias}_pychopper.tsv - workflow-glue generate_pychopper_stats --data ${meta.alias}_pychopper.tsv --output . + workflow-glue generate_pychopper_stats --data pychopper.tsv --output . - # Add sample id colum - sed "1s/\$/\tsample_id/; 1 ! s/\$/\t${meta.alias}/" ${meta.alias}_pychopper.tsv > tmp - mv tmp ${meta.alias}_pychopper.tsv + # Add sample id column + sed "1s/\$/\tsample_id/; 1 ! s/\$/\t${meta.alias}/" pychopper.tsv > tmp + mv tmp pychopper.tsv mkdir "${meta.alias}_pychopper_output/" @@ -274,9 +273,9 @@ process merge_gff_bundles{ input: tuple val(sample_id), path (gff_bundle) output: - tuple val(sample_id), path('*.gff'), emit: gff + tuple val(sample_id), path("${sample_id}.gff"), emit: gff script: - def merged_gff = "transcripts_${sample_id}.gff" + def merged_gff = "${sample_id}.gff" """ echo '##gff-version 2' >> $merged_gff; echo '#pipeline-nanopore-isoforms: stringtie' >> $merged_gff; @@ -303,13 +302,12 @@ process run_gffcompare{ tuple val(sample_id), path(query_annotation) path ref_annotation output: - tuple val(sample_id), path("${sample_id}_gffcompare"), - emit: gffcmp_dir + tuple val(sample_id), path("${sample_id}"), emit: gffcmp_dir path ("${sample_id}_annotated.gtf"), emit: gtf tuple val(sample_id), path("${sample_id}_transcripts_table.tsv"), emit: isoforms_table script: - def out_dir = "${sample_id}_gffcompare" + def out_dir = "${sample_id}" """ mkdir $out_dir echo "Doing comparison of reference annotation: ${ref_annotation} and the query annotation" @@ -320,14 +318,14 @@ process run_gffcompare{ workflow-glue generate_tracking_summary --tracking $out_dir/str_merged.tracking \ --output_dir ${out_dir} --annotation ${ref_annotation} - mv *.tmap $out_dir - mv *.refmap $out_dir - cp ${out_dir}/str_merged.annotated.gtf ${sample_id}_annotated.gtf + mv *.tmap "${out_dir}" + mv *.refmap "${out_dir}" + cp "${out_dir}/str_merged.annotated.gtf" "${sample_id}_annotated.gtf" # Make an isoform table for report and user output. workflow-glue make_isoform_table \ - --sample_id $sample_id \ - --gffcompare_dir "${sample_id}_gffcompare" + --sample_id "${sample_id}" \ + --gffcompare_dir "${out_dir}" """ } @@ -398,10 +396,9 @@ process makeReport { path "params.json" path "pychopper_report/*" path"jaffal_csv/*" - val sample_ids path "per_read_stats/?.gz" path "aln_stats/*" - path gffcmp_dir + path "gffcmp_dir/*" path "gff_annotation/*" path "de_report/*" path "seqkit/*" @@ -415,8 +412,6 @@ process makeReport { path ("filtered_transcript_counts_with_genes.tsv"), emit: filtered, optional: true path ("all_gene_counts.tsv"), emit: gene_counts, optional: true shell: - // Convert the sample_id arrayList. - sids = new BlankSeparatedList(sample_ids) report_name = "wf-transcriptomes-report.html" ''' if [ -f "de_report/OPTIONAL_FILE" ]; then @@ -430,10 +425,10 @@ process makeReport { else OPT_GFF_ANNOTATION="--gff_annotation gff_annotation/*" fi - if [ -f "OPTIONAL_FILE" ]; then + if [ -f "gffcmp_dir/OPTIONAL_FILE" ]; then OPT_GFFCMP_DIR="" else - OPT_GFFCMP_DIR="--gffcompare_dir !{gffcmp_dir}" + OPT_GFFCMP_DIR="--gffcompare_dir gffcmp_dir/" fi if [ -f "jaffal_csv/OPTIONAL_FILE" ]; then OPT_JAFFAL_CSV="" @@ -460,7 +455,6 @@ process makeReport { --params params.json \ ${OPT_ALN} \ ${OPT_PC_REPORT} \ - --sample_ids !{sids} \ --stats per_read_stats/* \ ${OPT_GFF_ANNOTATION} \ ${OPT_ISO_TABLE} \ @@ -596,7 +590,7 @@ workflow pipeline { } else{ full_len_reads = input_reads.map{ meta, reads -> [meta.alias, reads]} - pychopper_report = file("$projectDir/data/OPTIONAL_FILE") + pychopper_report = OPTIONAL_FILE } if (params.transcriptome_source != "precomputed"){ build_minimap_index(ref_genome) @@ -624,7 +618,7 @@ workflow pipeline { optional_channel = Channel.fromPath("$projectDir/data/OPTIONAL_FILE") gff_tuple = merge_gff_bundles.out.gff.combine(optional_channel) gff_compare = OPTIONAL_FILE - isoforms_table = file("$projectDir/data/OPTIONAL_FILE") + isoforms_table = OPTIONAL_FILE } // For reference based assembly, there is only one reference // So map this reference to all sample_ids @@ -638,10 +632,10 @@ workflow pipeline { results = results.concat(assembly.bam.map {sample_id, bam, bai -> [bam, bai]}.flatten()) } else{ - gff_compare = file("$projectDir/data/OPTIONAL_FILE") - isoforms_table = file("$projectDir/data/OPTIONAL_FILE") - merge_gff = file("$projectDir/data/OPTIONAL_FILE") - assembly_stats = file("$projectDir/data/OPTIONAL_FILE") + gff_compare = OPTIONAL_FILE + isoforms_table = OPTIONAL_FILE + merge_gff = OPTIONAL_FILE + assembly_stats = OPTIONAL_FILE use_ref_ann = false } @@ -649,7 +643,7 @@ workflow pipeline { gene_fusions(full_len_reads, jaffal_refBase, jaffal_genome, jaffal_annotation) jaffal_out = gene_fusions.out.results_csv.map{ alias, csv -> csv}.collectFile(keepHeader: true, name: 'jaffal.csv') } else{ - jaffal_out = file("$projectDir/data/OPTIONAL_FILE") + jaffal_out = OPTIONAL_FILE } if (params.de_analysis){ @@ -674,8 +668,8 @@ workflow pipeline { de_outputs = de.de_outputs counts = de.counts } else{ - de_report = file("$projectDir/data/OPTIONAL_FILE") - count_transcripts_file = file("$projectDir/data/OPTIONAL_FILE") + de_report = OPTIONAL_FILE + count_transcripts_file = OPTIONAL_FILE } makeReport( @@ -683,7 +677,6 @@ workflow pipeline { workflow_params, pychopper_report, jaffal_out, - input_reads.map{ meta, fastq -> meta.alias}.collect(), per_read_stats, assembly_stats, gff_compare, @@ -724,8 +717,8 @@ workflow pipeline { if (params.de_analysis){ de_results = report.concat( transcriptome, de_outputs.flatten(), counts.flatten(), - makeReport.out.results_dge, makeReport.out.tpm, - makeReport.out.filtered, makeReport.out.unfiltered, + makeReport.out.results_dge, makeReport.out.tpm, + makeReport.out.filtered, makeReport.out.unfiltered, makeReport.out.gene_counts) // Output de_analysis results in the dedicated directory. results = results.concat(de_results.map{ [it, "de_analysis"] }) @@ -768,7 +761,7 @@ workflow { error = "--ref_genome: File doesn't exist, check path." } }else { - ref_genome = file("$projectDir/data/OPTIONAL_FILE") + ref_genome = OPTIONAL_FILE } if (params.containsValue("denovo")) { error = "Denovo transcriptome source is no longer supported. Please use the reference-guided or precomputed options." @@ -782,7 +775,7 @@ workflow { } use_ref_ann = true }else{ - ref_annotation= file("$projectDir/data/OPTIONAL_FILE") + ref_annotation= OPTIONAL_FILE use_ref_ann = false } if (params.jaffal_refBase){ @@ -793,7 +786,7 @@ workflow { }else{ jaffal_refBase = null } - ref_transcriptome = file("$projectDir/data/OPTIONAL_FILE") + ref_transcriptome = OPTIONAL_FILE if (params.ref_transcriptome){ log.info("Reference Transcriptome provided will be used for differential expression.") ref_transcriptome = file(params.ref_transcriptome, type:"file") diff --git a/nextflow.config b/nextflow.config index f738b1c..338e853 100644 --- a/nextflow.config +++ b/nextflow.config @@ -116,7 +116,7 @@ manifest { description = 'Transcriptome analysis including gene fusions, differential expression as well as assembly and annotation of cDNA and direct RNA sequencing data.' mainScript = 'main.nf' nextflowVersion = '>=23.04.2' - version = 'v1.0.0' + version = 'v1.1.0' } epi2melabs { diff --git a/subworkflows/JAFFAL/gene_fusions.nf b/subworkflows/JAFFAL/gene_fusions.nf index 795106d..b186740 100644 --- a/subworkflows/JAFFAL/gene_fusions.nf +++ b/subworkflows/JAFFAL/gene_fusions.nf @@ -17,7 +17,7 @@ process jaffal{ # JAFFAL exists with status code 1 when there's 0 fusion hits. Prevent this with '||:' $params.jaffal_dir/tools/bin/bpipe run \ - -n $params.threads \ + -n "${task.cpus}" \ -p jaffa_output="\$JAFFAOUT/" \ -p refBase=$refBase \ -p genome=$genome \ diff --git a/subworkflows/reference_assembly.nf b/subworkflows/reference_assembly.nf index aa8f8f5..412759b 100644 --- a/subworkflows/reference_assembly.nf +++ b/subworkflows/reference_assembly.nf @@ -30,6 +30,11 @@ process map_reads{ | samtools sort --write-index -@ 1 -o "${sample_id}_reads_aln_sorted.bam##idx##${sample_id}_reads_aln_sorted.bam.bai" - ; ((cat "${sample_id}_reads_aln_sorted.bam" | seqkit bam -s -j 1 - 2>&1) | tee ${sample_id}_read_aln_stats.tsv ) || true + # Add sample id header and column + sed "s/\$/${sample_id}/" "${sample_id}_read_aln_stats.tsv" \ + | sed "1 s/${sample_id}/sample_id/" > tmp + mv tmp "${sample_id}_read_aln_stats.tsv" + if [[ -s "internal_priming_fail.tsv" ]]; then tail -n +2 "internal_priming_fail.tsv" | awk '{print ">" \$1 "\\n" \$4 }' - > "context_internal_priming_fail_start.fasta"