From c1f784ad04b5812c59bbcf5f4c935ed1b033d343 Mon Sep 17 00:00:00 2001 From: Sarah Griffiths Date: Mon, 31 Mar 2025 09:20:35 +0000 Subject: [PATCH] fix memory issue by replacing seq alignment stats with flagstat [CW-5858] --- CHANGELOG.md | 2 + bin/workflow_glue/de_plots.py | 69 +++++++++---------------- bin/workflow_glue/report.py | 4 +- main.nf | 6 +-- subworkflows/differential_expression.nf | 10 ++-- 5 files changed, 37 insertions(+), 54 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5c17cdc..5013138 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,8 +12,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - A common user issue is providing a ref_annotation and ref_genome parameter that have mismatched reference IDs, which causes the DE_analysis to fail. The workflow will now do an upfront check and give an error message if no overlap is found or a warning if some IDs are present in one file but not in the other. - Reconciled workflow with wf-template v5.5.0. - Sort the columns and rows of the gene and transcript count files. +- DE_analysis alignment summary stats table no longer includes MAPQ or quality scores. MAPQ is not relevant for transcript alignment and quality scores are already available in the read summary section of the report. ### Fixed - `all_gene_counts.tsv` contained the DE counts results. +- Reduced memory usage of the report workflow process. - The merged transcriptome generated for differential expression analysis now only contains the exons and not the full genomic sequence. - Output the gene name annotated differential expression analysis count files only. - Only use full length reads in the differential expression analysis. diff --git a/bin/workflow_glue/de_plots.py b/bin/workflow_glue/de_plots.py index 2a313c8..bc41845 100644 --- a/bin/workflow_glue/de_plots.py +++ b/bin/workflow_glue/de_plots.py @@ -1,5 +1,6 @@ #!/usr/bin/env python """Create de report section.""" +import json import os from dominate.tags import h5, p @@ -7,48 +8,32 @@ from dominate.util import raw from ezcharts import scatterplot from ezcharts.components.ezchart import EZChart from ezcharts.layout.snippets import DataTable +from natsort import natsorted import numpy as np import pandas as pd -def parse_seqkit(fname): - """Get seqkit columns.""" - cols = { - 'Read': str, 'Ref': str, 'MapQual': int, 'Acc': float, 'ReadLen': int, - 'ReadAln': int, 'ReadCov': float, 'MeanQual': float, - 'IsSec': bool, 'IsSup': bool} - df = pd.read_csv(fname, sep="\t", dtype=cols, usecols=cols.keys()) - df['Clipped'] = df['ReadLen'] - df['ReadAln'] - df['Type'] = 'Primary' - df.loc[df['IsSec'], 'Type'] = 'Secondary' - df.loc[df['IsSup'], 'Type'] = 'Supplementary' - df["fname"] = os.path.basename(fname).rstrip(".seqkit.stats") - return df - - -def number_of_alignments(df, field_name): - """Group alignments for summary table.""" - grouped = df.groupby('fname').agg(**{ - field_name: ('Read', 'size'), - }) - return grouped.transpose() - - -def create_summary_table(df): - """Create summary table.""" - all_aln = number_of_alignments(df, "Read mappings") - primary = number_of_alignments(df.loc[df['Type'] == 'Primary'], "Primary") - secondary = number_of_alignments( - df.loc[df['Type'] == 'Secondary'], "Secondary") - supplementary = number_of_alignments( - df.loc[df['Type'] == 'Supplementary'], "Supplementary") - avg_acc = df.loc[df['Type'] == 'Primary'].groupby( - 'fname').agg(**{"Median Qscore": ('MeanQual', 'median'), }).transpose() - avg_mapq = df.loc[df['Type'] == 'Primary'].groupby( - 'fname').agg(**{"Median MAPQ": ('MapQual', 'median'), }).transpose() - return pd.concat([ - all_aln, primary, secondary, supplementary, - avg_acc, avg_mapq]) +def flagstats_df(flagstats_reports): + """Flag stats alignment dataframe.""" + flagstats_dic = {} + for flagstat in flagstats_reports.iterdir(): + with open(flagstat, "r") as f: + data = json.load(f) + data = data["QC-passed reads"] + flagstats = [ + 'mapped', 'primary mapped', 'secondary', 'supplementary'] + per_sample_flagstats = {key: data.get(key) for key in flagstats} + sample = os.path.basename(flagstat).split(".")[0] + flagstats_dic[sample] = per_sample_flagstats + alignment_summary_df = pd.DataFrame(flagstats_dic) + alignment_summary_df = alignment_summary_df[ + natsorted(alignment_summary_df.columns)] + alignment_summary_df.index = [ + "Total Read Mappings", + "Primary", "Secondary", + "Supplementary"] + alignment_summary_df.index.name = "Statistic" + return alignment_summary_df def dexseq_section(dexseq_file, tr_id_to_gene_name, tr_id_to_gene_id, pval_thresh): @@ -239,7 +224,7 @@ def get_translations(gtf): def de_section( annotation, dge, dexseq, dtu, tpm, report, filtered, unfiltered, - gene_counts, aln_stats_dir, pval_threshold=0.01): + gene_counts, flagstats_dir, pval_threshold=0.01): """Differential expression sections.""" with (report.add_section("Differential expression", "DE")): @@ -255,13 +240,9 @@ def de_section( Find the full sequences of any transcripts in the final_non_redundant_transcriptome.fasta file. """) - alignment_stats = pd.concat([parse_seqkit(f) for f in aln_stats_dir.iterdir()]) - alignment_summary_df = create_summary_table(alignment_stats) - alignment_summary_df = alignment_summary_df.fillna(0).applymap(np.int64) + alignment_summary_df = flagstats_df(flagstats_dir) h5("Alignment summary stats") - alignment_summary_df.index.name = "statistic" DataTable.from_pandas(alignment_summary_df, use_index=True) - salmon_table(tpm) # Get translations for adding gene names to tables diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py index 8ad510e..14692f4 100644 --- a/bin/workflow_glue/report.py +++ b/bin/workflow_glue/report.py @@ -314,7 +314,7 @@ def transcriptome_summary(report, summaries_dir): df_table, use_index=False, searchable=False, paging=False) -def de_section(report, de_report_dir, de_aln_stats_dir, pval_threshold): +def de_section(report, de_report_dir, flagstats_dir, pval_threshold): """Make differential transcript expression section.""" dexseq = de_report_dir / "results_dexseq.tsv" dge = de_report_dir / "results_dge.tsv" @@ -337,7 +337,7 @@ def de_section(report, de_report_dir, de_aln_stats_dir, pval_threshold): filtered=filtered, unfiltered=unfiltered, gene_counts=gene_counts, - aln_stats_dir=de_aln_stats_dir, + flagstats_dir=flagstats_dir, pval_threshold=pval_threshold ) diff --git a/main.nf b/main.nf index 5a26c71..dee81fe 100644 --- a/main.nf +++ b/main.nf @@ -757,10 +757,10 @@ workflow pipeline { de = differential_expression(transcriptome, full_len_reads.map{ sample_id, reads -> [[alias:sample_id], reads]}, sample_sheet, gtf) de_report = de.all_de de_outputs = de.de_outputs - count_transcripts_file = de.count_transcripts + de_alignment_stats = de.de_alignment_stats } else{ de_report = OPTIONAL_FILE - count_transcripts_file = OPTIONAL_FILE + de_alignment_stats = OPTIONAL_FILE } // get metadata and stats files, keeping them ordered (could do with transpose I suppose) @@ -777,7 +777,7 @@ workflow pipeline { software_versions, workflow.manifest.version, workflow_params, - count_transcripts_file, + de_alignment_stats, pychopper_report, assembly_stats, gff_compare, diff --git a/subworkflows/differential_expression.nf b/subworkflows/differential_expression.nf index 86e6370..5e1c97b 100644 --- a/subworkflows/differential_expression.nf +++ b/subworkflows/differential_expression.nf @@ -16,16 +16,14 @@ process count_transcripts { // library type is specified as forward stranded (-l SF) as it should have either been through pychopper or come from direct RNA reads. label "isoforms" cpus params.threads - memory "16 GB" + memory "31 GB" input: tuple val(meta), path(bam), path(ref_transcriptome) output: path "*transcript_counts.tsv", emit: counts - path "*seqkit.stats", emit: seqkit_stats """ salmon quant --noErrorModel -p "${task.cpus}" -t "${ref_transcriptome}" -l SF -a "${bam}" -o counts mv counts/quant.sf "${meta.alias}.transcript_counts.tsv" - seqkit bam "${bam}" 2> "${meta.alias}.seqkit.stats" """ } @@ -147,10 +145,12 @@ process map_transcriptome{ tuple val(meta), path (fastq_reads), path(index) output: tuple val(meta), path("${meta.alias}_reads_aln_sorted.bam"), emit: bam + path("${meta.alias}.flagstat.stats"), emit: align_stats """ minimap2 -t ${task.cpus} -ax splice -uf -p 1.0 "${index}" "${fastq_reads}" \ | samtools view -Sb > "output.bam" samtools sort -@ ${task.cpus} "output.bam" -o "${meta.alias}_reads_aln_sorted.bam" + samtools flagstat -O json "${meta.alias}_reads_aln_sorted.bam" > "${meta.alias}.flagstat.stats" """ } @@ -179,9 +179,9 @@ workflow differential_expression { // Concat files required to be output to user without any changes de_outputs_concat = analysis.cpm.concat(plotResults.out.dtu_plots, analysis.dge_pdf, analysis.dge_tsv, analysis.dtu_gene, analysis.dtu_transcript, analysis.dtu_stageR, analysis.dtu_pdf, merged_TPM).collect() - count_transcripts_file = count_transcripts.out.seqkit_stats.collect() + collected_de_alignment_stats = mapped.align_stats.collect() emit: all_de = de_report - count_transcripts = count_transcripts_file + de_alignment_stats = collected_de_alignment_stats de_outputs = de_outputs_concat }