diff --git a/CHANGELOG.md b/CHANGELOG.md index e4ff3b8..b79e883 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,20 +4,22 @@ 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.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -# Unreleased -## Added -- --direct_rna option. +## [v0.1.2] +### Added +- direct_rna option +- Some extra error handling +- Minor report display improvements ## [v0.1.1] -## Fixed +### Fixed - Incorrect numbers and of transcripts caused by merging gff files with same gene and transcript ids - Error handling in de novo pipeline. Skip clusters in build_backbones that cause an isONclust2 error - Several small fixes in report plotting ## [v0.1.0] -## Added +### Added - Added the denovo pipeline -## Changed +### Changed - Updates to the report plots ## [v0.0.1] diff --git a/bin/report.py b/bin/report.py index 8d37efb..bbc2024 100755 --- a/bin/report.py +++ b/bin/report.py @@ -7,12 +7,12 @@ import math from pathlib import Path from aplanat import bars, hist, lines -from aplanat.components import fastcat from aplanat.components import simple as scomponents +from aplanat.components.fastcat import read_length_plot, read_quality_plot from aplanat.report import WFReport from aplanat.util import Colors from bokeh.layouts import gridplot -from bokeh.models import ColumnDataSource, Panel, Tabs +from bokeh.models import ColumnDataSource, Legend, Panel, Tabs from bokeh.models.widgets import DataTable, TableColumn from bokeh.palettes import Category10_10 from bokeh.plotting import figure @@ -253,6 +253,8 @@ def grouped_bar(df, title="", tilted_xlabs=False): dodge_increment = abs(dodge_range[0] - dodge_range[1]) / \ (len(df.columns) - 1) + legend_it = [] + if tilted_xlabs: p.xaxis.major_label_orientation = math.pi / 4 @@ -270,15 +272,18 @@ def grouped_bar(df, title="", tilted_xlabs=False): i += 1 width = df.size / 60 - p.vbar( + v = p.vbar( x=dodge('x_groups', current_dodge, range=p.x_range), top=col, - width=width, source=source, color=color, legend_label=col) + width=width, source=source, color=color) current_dodge += dodge_increment + legend_it.append([col, [v]]) + legend = Legend(items=legend_it) + p.add_layout(legend, 'right') p.x_range.range_padding = 0.1 p.xgrid.grid_line_color = None p.legend.location = "top_left" - p.legend.orientation = "horizontal" + p.legend.orientation = "vertical" return p @@ -319,6 +324,7 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids): ''') tabs = [] + gff_fails = False for id_, dir_ in zip(sample_ids, gffcompare_outdirs): stats, _, miss, novel, total = \ parse_gffcmp_stats(dir_ / 'str_merged.stats') @@ -333,16 +339,17 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids): child=gridplot( [bar_totals, bar_performance, bar_missed, bar_novel], ncols=2, width=350, height=260), title=id_)) - - cover_panel = Tabs(tabs=tabs) - section.markdown(gffcompare_md) - section.plot(cover_panel) else: - gffcompare_md = (''' - __Warning__: Gffcompare summary cannot be shown. - This could be due to incompatible reference fasta and gff - files.''') - section.markdown(gffcompare_md) + gff_fails = True + + if gff_fails: + gffcompare_md += (''' + __Warning__: Some gffcompare summary cannot be shown. + This could be due to incompatible reference fasta and gff + files.''') + cover_panel = Tabs(tabs=tabs) + section.markdown(gffcompare_md) + section.plot(cover_panel) names = { '=': 'ExactMatch:=', @@ -598,9 +605,9 @@ def transcript_table(report, df_tmaps, covr_threshold): ### Query transcript table Low coverage transcripts are removed to speed up the table viewing.
- This can be set with the parameter `transcript_table_cov_thresh` in the - config. - ''') + Coverage threshold can be set with the parameter + `transcript_table_cov_thresh`. + ''') df = df_tmaps.drop( columns=[ @@ -614,6 +621,12 @@ def transcript_table(report, df_tmaps, covr_threshold): df.sort_values('cov', ascending=True, inplace=True) counts = list(range(len(df))) + # Filter on coverage threshold + df = df[df['cov'] >= covr_threshold] + if len(df) < 200: # Min size of table should be 200. Don't filter + df = df.sort_values('cov', ascending=False).iloc[:, 0:200] + covr_threshold = 0 + # Keep Isoforms with coverage > threshold vline_x = np.argmax(df['cov'] > covr_threshold) vline_y = [0, df['cov'].max()] @@ -629,10 +642,6 @@ def transcript_table(report, df_tmaps, covr_threshold): colors=['blue', 'red']) section.plot(cov_plt) - # Filter on coverage threshold - df = df[df['cov'] >= covr_threshold] - if len(df) < 200: # Min size of table should be 200 - df = df.sort_values('cov', ascending=False).iloc[:, 0:200] # Make a column of number of isoforms in parent gene gb = df.groupby(['ref_gene_id', 'sample_id']).count() @@ -707,6 +716,8 @@ def transcriptome_summary(report, gffs, sample_ids, denovo=False): bar_isos = hist.histogram( [isoforms_per_gene], colors=[Colors.cerulean], title="isoforms per gene") + bar_isos.xaxis.axis_label = "Num. isoforms" + bar_isos.yaxis.axis_label = "Num. genes" bar_isos.xaxis.major_label_orientation = math.pi / 2.8 plots.append(bar_isos) @@ -723,6 +734,8 @@ def transcriptome_summary(report, gffs, sample_ids, denovo=False): fig = figure(title="Exons per transcript") fig.vbar( x, top=list(y), color=Colors.cerulean) + fig.xaxis.axis_label = 'Num. exons' + fig.yaxis.axis_label = 'Num. genes' fig.xaxis.major_label_orientation = math.pi / 2.8 plots.append(fig) @@ -765,6 +778,23 @@ def load_sample_data(files, sample_ids, read_func=None): return df_ +def seq_stats_tabs(report, sample_ids, stats): + """Make tabs of sequence summaries by sample.""" + tabs = [] + for id_, summ in sorted(zip(sample_ids, stats)): + df_sum = pd.read_csv(summ, index_col=False, sep='\t') + rlp = read_length_plot(df_sum) + rqp = read_quality_plot(df_sum) + grid = gridplot( + [rlp, rqp], ncols=2, sizing_mode="stretch_width") + + tabs.append(Panel(child=grid, title=id_)) + section = report.add_section() + section.markdown(""" + ### Sequence summaries""") + section.plot(Tabs(tabs=tabs)) + + def main(): """Run the entry point.""" parser = argparse.ArgumentParser() @@ -815,12 +845,7 @@ def main(): revision=args.revision, commit=args.commit) # Add reads summary section - for id_, summ in zip(sample_ids, args.summaries): - report.add_section( - section=fastcat.full_report( - [summ], - header='#### Read stats: {}'.format(id_) - )) + seq_stats_tabs(report, args.sample_ids, args.summaries) if args.alignment_stats is not None: df_aln_stats = load_sample_data(args.alignment_stats, sample_ids) diff --git a/main.nf b/main.nf index af60b1c..2f230fd 100644 --- a/main.nf +++ b/main.nf @@ -130,21 +130,32 @@ process split_bam{ output: tuple val(sample_id), path('*.bam'), emit: bundles script: - if (params["bundle_min_reads"] != false) - """ + """ + n=`samtools view -c $bam` + if [[ n -lt 1 ]] + then + echo 'There are no reads mapping for $sample_id. Exiting!' + exit 1 + fi + + re='^[0-9]+\$' + + if [[ $params.bundle_min_reads =~ \$re ]] + then + echo "Bundling up the bams" seqkit bam -j ${params.threads} -N ${params.bundle_min_reads} ${bam} -o bam_bundles/ - mv bam_bundles/* . let i=1 - for b in *.bam; do - newname="${sample_id}_batch_\${i}.bam" - mv \$b \$newname - ((i++)) + for b in bam_bundles/*.bam; do + echo \$b + newname="${sample_id}_batch_\${i}.bam" + mv \$b \$newname + ((i++)) done - """ else - """ + echo 'no bundling' ln -s ${bam} ${sample_id}_batch_1.bam - """ + fi + """ } diff --git a/nextflow.config b/nextflow.config index 33f4bf6..f4b52c6 100644 --- a/nextflow.config +++ b/nextflow.config @@ -23,7 +23,7 @@ params { sample = null sample_sheet = null sanitize_fastq = false - wfversion = "v0.1.1" + wfversion = "v0.1.2" aws_image_prefix = null aws_queue = null report_name = "report" @@ -34,7 +34,7 @@ params { schema_ignore_params = 'show_hidden_params,validate_params,monochrome_logs,aws_queue,aws_image_prefix,wfversion' // Process cDNA reads using pychopper, turn off for direct RNA: - direct_rna = true + direct_rna = false // Options passed to pychopper: pychopper_opts = "-m edlib" diff --git a/nextflow_schema.json b/nextflow_schema.json index 7e194eb..4e14e9f 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -280,7 +280,7 @@ }, "wfversion": { "type": "string", - "default": "v0.1.1", + "default": "v0.1.2", "hidden": true }, "monochrome_logs": {