Merge branch 'CW-5858' into 'dev'

fix memory issue by replacing seq alignment stats with flagstat [CW-5858]

Closes CW-5858

See merge request epi2melabs/workflows/wf-transcriptomes!208

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.

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

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
     )
 

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,

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
 }