diff --git a/bin/report.py b/bin/report.py index 227eb1c..a15431d 100755 --- a/bin/report.py +++ b/bin/report.py @@ -3,7 +3,6 @@ import argparse from collections import Counter, defaultdict, OrderedDict -from functools import reduce import math from pathlib import Path @@ -19,7 +18,6 @@ from bokeh.palettes import Category10_10 from bokeh.plotting import figure from bokeh.transform import dodge import gffutils -from jinja2 import Template import numpy as np import pandas as pd import sigfig @@ -50,30 +48,6 @@ def _hbar(y, right, title='', fig_height=300, fig_width=300, return fig -class Table: - """A table report component. - - Adapted from aplanat - """ - - def __init__(self, template, data_frame, index, table_id, **kwargs): - """Initialize table component. - - :param dataframe: dataframe to turn in to simple table. - """ - template = Template(template) - - for key, val in kwargs.items(): - if isinstance(val, bool): - kwargs[key] = str(val).lower() - - self.div = template.render(dataframe=data_frame.to_html( - table_id=table_id, - index=index), - table_id=table_id, - kwargs=kwargs) - - def simple_hbar(df, y, right, title="", color=Colors.cerulean, fig_kwargs={}, plot_kwargs={}): """Create a simple barplot. @@ -461,33 +435,22 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids): tracking['Overlaps'].values.tolist(), tracking['Percent'].values.tolist(), title="{}".format(id_)) - # Edit for creating unified table - tracking.rename(columns={'sample_id': id_ + ' count'}, inplace=True) - tracking.drop(columns=['Count'], inplace=True) + tracking.drop(columns=['sample_id'], inplace=True) tracking_dfs.append(tracking) - # section.plot(grid) - - # Tracking table - df_class_table = reduce( - lambda left, right: pd.merge(left, right), tracking_dfs) - - desc = pd.Series(df_class_table.Overlaps.apply( + tracking['description'] = pd.Series(tracking.Overlaps.apply( lambda x: x.split(':')[0])) - df_class_table.insert(0, 'description', desc) - code = pd.Series(df_class_table.Overlaps.apply( + tracking['code'] = pd.Series(tracking.Overlaps.apply( lambda x: x.split(':')[1])) - df_class_table.insert(0, 'code', code) cols = [TableColumn(field=Ci, title=Ci, width=100) - for Ci in df_class_table.columns] + for Ci in tracking.columns] + track_table = DataTable(columns=cols, - source=ColumnDataSource(df_class_table), + source=ColumnDataSource(tracking), index_position=None, width=500) - - df_class_table.drop(columns=['Overlaps'], inplace=True) tabs.append(Panel( child=gridplot([track_bar, track_table], ncols=2), title=id_) ) @@ -654,7 +617,7 @@ def cluster_quality(cluster_qc_dir, report, sample_ids): section.plot(cover_panel) -def transcript_table(report, df_tmaps, covr_threshold, table_template): +def transcript_table(report, df_tmaps, covr_threshold): """Create searchable table of transcripts.""" section = report.add_section() @@ -701,11 +664,7 @@ def transcript_table(report, df_tmaps, covr_threshold, table_template): df.sort_values('parent gene iso num', inplace=True, ascending=True) - with open(table_template, 'r') as fh: - tabletempl = fh.read() - - bigtable = Table(tabletempl, df, index=False, table_id='bigtable') - section._add_item(bigtable.div) + section.table(df, index=False) def tanscriptome_summary(report, gffs, sample_ids, denovo=False): @@ -868,12 +827,6 @@ def main(): parser.add_argument( "--sample_ids", required=True, nargs='+', help="List of sample ids") - parser.add_argument( - "--report_template", required=True, - help="Jinja template") - parser.add_argument( - "--table_template", required=True, - help="Template for big transcript table") parser.add_argument( "--transcript_table_cov_thresh", required=False, type=int, default=50, help="Isoforms without this support will be excluded from the table") @@ -932,14 +885,8 @@ def main(): if len(pc_df) > 0: pychopper_plots(report, pc_df) - with open(args.report_template, "r") as fh: - reptempl = fh.read() - - report.template = Template(reptempl) - if df_tmaps is not None: - transcript_table(report, df_tmaps, args.transcript_table_cov_thresh, - args.table_template) + transcript_table(report, df_tmaps, args.transcript_table_cov_thresh) if args.cluster_qc_dirs is not None: cluster_quality(args.cluster_qc_dirs, report, sample_ids) diff --git a/bin/report_template.html b/bin/report_template.html deleted file mode 100755 index 57fe286..0000000 --- a/bin/report_template.html +++ /dev/null @@ -1,42 +0,0 @@ - - - - - - - - wf-isoforms report - - - - - - - - - {{ resources }} - {{ script }} - - {{ bigtable_js }} - - -
-

{{ title }}

-

{{ lead }} - {{ div }} -

- {{ big_table }} -
- - diff --git a/bin/table_template.html b/bin/table_template.html deleted file mode 100755 index 1faf6d9..0000000 --- a/bin/table_template.html +++ /dev/null @@ -1,106 +0,0 @@ - - - - - -
Table loading
-{{ dataframe }} - - - diff --git a/evaluation/tests.sh b/evaluation/tests.sh old mode 100644 new mode 100755 index 515016f..59daaff --- a/evaluation/tests.sh +++ b/evaluation/tests.sh @@ -26,7 +26,6 @@ multisampledir="test_data/demultiplexed_fastq" #"--minimap2_opts '-uf --splice-flank=no'" results=() - OUTPUT=$1/denovo_multi_sample_no_ref_genome; nextflow run . --fastq $multisampledir $config --denovo --ref_genome test_data/SIRV_150601a.fasta -profile local --out_dir ${OUTPUT} -w ${OUTPUT}/workspace \ --sample_sheet test_data/sample_sheet -resume; diff --git a/main.nf b/main.nf index 6cd0250..e0d2329 100644 --- a/main.nf +++ b/main.nf @@ -42,7 +42,6 @@ process getVersions { python -c "import pysam; print(f'pysam,{pysam.__version__}')" >> versions.txt python -c "import aplanat; print(f'aplanat,{aplanat.__version__}')" >> versions.txt python -c "import pandas; print(f'pandas,{pandas.__version__}')" >> versions.txt - python -c "import sklearn; print(f'scikit-learn,{sklearn.__version__}')" >> versions.txt fastcat --version | sed 's/^/fastcat,/' >> versions.txt minimap2 --version | sed 's/^/minimap2,/' >> versions.txt samtools --version | head -n 1 | sed 's/ /,/' >> versions.txt @@ -87,15 +86,9 @@ process preprocess_reads { tuple val(sample_id), path('*.tsv'), emit: report script: """ - if [[ ${params.use_pychopper} == true ]]; - then - cdna_classifier.py -t ${params.threads} ${params.pychopper_opts} ${input_reads} ${sample_id}_full_length_reads.fq - mv cdna_classifier_report.tsv ${sample_id}_cdna_classifier_report.tsv - generate_pychopper_stats.py --data ${sample_id}_cdna_classifier_report.tsv --output . - else - ln -s `realpath $input_reads` "${sample_id}_full_length_reads.fq" - touch $sample_id}_cdna_classifier_report.tsv - fi + cdna_classifier.py -t ${params.threads} ${params.pychopper_opts} ${input_reads} ${sample_id}_full_length_reads.fq + mv cdna_classifier_report.tsv ${sample_id}_cdna_classifier_report.tsv + generate_pychopper_stats.py --data ${sample_id}_cdna_classifier_report.tsv --output . """ } @@ -267,8 +260,6 @@ process makeReport { label "isoforms" input: - path report_template - path table_template path versions path "params.json" val denovo @@ -289,8 +280,6 @@ process makeReport { def OPT_DENOVO = denovo ? "--denovo" : '' """ report.py --report $report_name \ - --report_template $report_template \ - --table_template $table_template \ --versions $versions \ --params params.json \ $OPT_ALN \ @@ -301,8 +290,6 @@ process makeReport { --gff_annotation $gff_annotation \ --transcript_table_cov_thresh $params.transcript_table_cov_thresh \ $OPT_DENOVO - - """ } @@ -311,43 +298,24 @@ process makeReport { // decoupling the publish from the process steps. process output { // publish inputs to output directory - label "isoforms" - publishDir "${params.results_dir}/${sample_id}", mode: 'copy', pattern: "*" - + publishDir "${params.out_dir}", mode: 'copy', pattern: "*" input: - tuple val(sample_id), path(fname) + path fname output: path fname """ echo "Writing output files" - echo $fname """ } -process output_report { - publishDir "${params.results_dir}", mode: 'copy', pattern: "*report.html" - - input: - path fname - output: - path fname - """ - echo "Copying report" - """ -} - - // workflow module workflow pipeline { take: reads ref_genome ref_annotation - report_template - table_template main: - - map_sample_ids_cls = {it -> + map_sample_ids_cls = {it -> /* Harmonize tuples output: tuple val(sample_id), path('*.gff') @@ -405,9 +373,7 @@ workflow pipeline { seq_for_transcriptome_build = sample_ids.flatten().combine(Channel.fromPath(params.ref_genome)) } - makeReport(report_template, - table_template, - software_versions, + makeReport(software_versions, workflow_params, params.denovo, summariseConcatReads.out.summary @@ -424,18 +390,22 @@ workflow pipeline { .join(run_gffcompare.out.gffcmp_dir) .join(seq_for_transcriptome_build)) - if (use_ref_ann){ + if (use_ref_ann){ results = preprocess_reads.out.report .concat(run_gffcompare.output.gffcmp_dir, m.stats, - get_transcriptome.out.flatMap(map_sample_ids_cls), - ) - } + get_transcriptome.out.flatMap(map_sample_ids_cls)) + .map {it -> it[1]} + .concat(makeReport.out.report) + + } if (!use_ref_ann && !params.denovo){ results = preprocess_reads.out.report .concat(m.stats, - get_transcriptome.out.flatMap(map_sample_ids_cls), - ) + get_transcriptome.out.flatMap(map_sample_ids_cls)) + .map {it -> it[1]} + .concat(makeReport.out.report) + } if (params.denovo){ results = m.cds @@ -443,18 +413,19 @@ workflow pipeline { seq_for_transcriptome_build, get_transcriptome.out.flatMap(map_sample_ids_cls), merge_gff_bundles.out.gff, - m.opt_qual_ch.flatMap {it -> + m.opt_qual_ch.flatMap {it -> l = [] for (x in it[1..-1]){ l.add(tuple(it[0], x)) } return l }) + .map {it -> it[1]} + .concat(makeReport.out.report) } emit: results - report telemetry = workflow_params } @@ -463,10 +434,6 @@ WorkflowMain.initialise(workflow, params, log) workflow { start_ping() - params.results_dir = "${params.out_dir}/output" - - report_template = file("$projectDir/bin/report_template.html") - table_template = file("$projectDir/bin/table_template.html") fastq = file(params.fastq, type: "file") @@ -510,10 +477,9 @@ workflow { params.fastq, params.out_dir, params.sample, params.sample_sheet, params.sanitize_fastq ) - pipeline(reads, ref_genome, ref_annotation, report_template, table_template) + pipeline(reads, ref_genome, ref_annotation) output(pipeline.out.results) - output_report(pipeline.out.report) end_ping(pipeline.out.telemetry) } diff --git a/nextflow.config b/nextflow.config index afa7082..7962535 100644 --- a/nextflow.config +++ b/nextflow.config @@ -15,17 +15,15 @@ params { fastq = null ref_genome = false ref_annotation = null - // Process cDNA reads using pychopper, turn off for direct RNA: - use_pychopper = true threads = 4 // Thresholds for viewing isoforms in report table transcript_table_cov_thresh = 50 - out_dir = null + out_dir = "output" sample = null sample_sheet = null sanitize_fastq = false - wfversion = "v0.1.0" + wfversion = "v0.1.1" aws_image_prefix = null aws_queue = null report_name = "report" diff --git a/nextflow_schema.json b/nextflow_schema.json index 5ea1c32..392c5c8 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -61,11 +61,6 @@ "type": "integer", "default": 4 }, - "use_pychopper": { - "type": "boolean", - "description": "Use pychopper to preprcess reads", - "default": true - }, "pychopper_opts": { "type": "string", "description": "Extra pychopper opts", diff --git a/reference_assembly.nf b/reference_assembly.nf index e0da46a..d88767b 100644 --- a/reference_assembly.nf +++ b/reference_assembly.nf @@ -24,8 +24,8 @@ process map_reads{ minimap2 -t ${params.threads} -ax splice ${params.minimap2_opts} ${index} "reads.fa"\ | samtools view -q ${params.minimum_mapping_quality} -F 2304 -Sb -\ | seqkit bam -j ${params.threads} -x -T '${ContextFilter}' -\ - | samtools sort -@ ${params.threads} -o "${sample_id}_reads_aln_sorted.bam" - \ - | ((seqkit bam -s -j ${params.threads} - 2>&1) | tee ${sample_id}_read_aln_stats.tsv ) || true + | samtools sort -@ ${params.threads} -o "${sample_id}_reads_aln_sorted.bam" - ; + ((cat "${sample_id}_reads_aln_sorted.bam" | seqkit bam -s -j ${params.threads} - 2>&1) | tee ${sample_id}_read_aln_stats.tsv ) || true if [[ -s "internal_priming_fail.tsv" ]]; then