From d477a747e6e1c93809e3fd93ed7121519c5af1e7 Mon Sep 17 00:00:00 2001 From: Neil Horner Date: Mon, 8 Aug 2022 13:48:15 +0000 Subject: [PATCH] Switch isoform table thresholding to just displaying to n rows based on coverage --- bin/report.py | 63 +++++++++++++++++--------------------------- main.nf | 2 +- nextflow.config | 2 +- nextflow_schema.json | 6 ++--- 4 files changed, 29 insertions(+), 44 deletions(-) diff --git a/bin/report.py b/bin/report.py index 47dd946..b31937c 100755 --- a/bin/report.py +++ b/bin/report.py @@ -6,7 +6,7 @@ from collections import Counter, defaultdict, OrderedDict import math from pathlib import Path -from aplanat import bars, hist, lines +from aplanat import bars, hist from aplanat.components import simple as scomponents from aplanat.components.fastcat import read_length_plot, read_quality_plot from aplanat.report import WFReport @@ -594,21 +594,25 @@ def cluster_quality(cluster_qc_dir, report, sample_ids): section.plot(cover_panel) -def transcript_table(report, df_tmaps, covr_threshold): - """Create searchable table of transcripts.""" +def transcript_table(report, df_tmaps, max_rows): + """Create searchable table of transcripts. + + :param df_tmaps: pd.DataFrame of all gffcomapre `.tmap` files from all + samples + """ section = report.add_section() - # Should we put data from each sample into it's own table or have it + # Should wße put data from each sample into it's own table or have it # all in single table and sample_id column? Currently it's the latter # drop some columns for the big table and do some filtering section.markdown(''' ### Isoforms table - Low coverage transcripts are removed to speed up the table viewing.
- Coverage threshold can be set with the parameter - `transcript_table_cov_thresh`. - ''') + Table interactivity can be slow if too many isoforms are loaded.
+ The number of isoform rows to load in this table can be set with + `isoform_table_nrows`. It is currently set to `{}` + '''.format(max_rows)) df = df_tmaps.drop( columns=[ @@ -619,31 +623,6 @@ def transcript_table(report, df_tmaps, covr_threshold): section.markdown("No transcripts found") return - 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()] - - cov_plt = lines.line( - [counts, [vline_x, vline_x]], # x-values - [df['cov'].values.tolist(), vline_y], # y-values - title=( - "Read Coverage. Threshold = {}x coverage".format( - covr_threshold) - ), x_axis_label='Num Isoforms', - y_axis_label='Coverage', - colors=['blue', 'red']) - - section.plot(cov_plt) - # Make a column of number of isoforms in parent gene gb = df.groupby(['ref_gene_id', 'sample_id']).count() # gb = gb.set_index(['ref_gene_id', 'sample_id']) @@ -654,7 +633,13 @@ def transcript_table(report, df_tmaps, covr_threshold): # Uncalssified transcritps should not be lumped togetehr df.loc[df.class_code == 'u', 'parent gene iso num'] = None - df.sort_values('parent gene iso num', inplace=True, ascending=True) + # Keep top n rows with most coverage + df.sort_values('cov', ascending=False, inplace=True) + df['cov'] = df['cov'].astype(int) + df = df.iloc[0: max_rows, :] + + # Sort by transcripts with highest isoform diversity + df.sort_values('parent gene iso num', inplace=True, ascending=False) section.table(df, index=False) @@ -664,7 +649,7 @@ def transcriptome_summary(report, gffs, sample_ids, denovo=False): Plot transcriptome summaries. Some of this data is available via gffcompare output, but the de novo - pipeline skips that, so we do it al here. + pipeline skips that, so we do it all here. We do not report exon number for the denovo assembly yet. This is because in this case, the gff annotation is generated by aligning to the CDS not @@ -718,7 +703,7 @@ def transcriptome_summary(report, gffs, sample_ids, denovo=False): bar_isos = hist.histogram( [isoforms_per_gene], colors=[Colors.cerulean], - title="isoforms per gene") + title="isoforms per gene", binwidth=1) bar_isos.xaxis.axis_label = "Num. isoforms" bar_isos.yaxis.axis_label = "Num. genes" @@ -864,8 +849,8 @@ def main(): "--sample_ids", required=True, nargs='+', help="List of sample ids") parser.add_argument( - "--transcript_table_cov_thresh", required=False, type=int, default=50, - help="Isoforms without this support will be excluded from the table") + "--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") @@ -912,7 +897,7 @@ def main(): report.write(args.report) if df_tmaps is not None: - transcript_table(report, df_tmaps, args.transcript_table_cov_thresh) + transcript_table(report, df_tmaps, args.isoform_table_nrows) if args.cluster_qc_dirs is not None: cluster_quality(args.cluster_qc_dirs, report, sample_ids) diff --git a/main.nf b/main.nf index eea88ba..d868e54 100644 --- a/main.nf +++ b/main.nf @@ -309,7 +309,7 @@ process makeReport { --summaries $seq_summaries \ --gffcompare_dir $gffcmp_dir \ --gff_annotation $gff_annotation \ - --transcript_table_cov_thresh $params.transcript_table_cov_thresh \ + --isoform_table_nrows $params.isoform_table_nrows \ $OPT_JAFFAL_CSV \ $OPT_DENOVO """ diff --git a/nextflow.config b/nextflow.config index fa1aa20..e058a4c 100644 --- a/nextflow.config +++ b/nextflow.config @@ -17,7 +17,7 @@ params { ref_annotation = null threads = 4 // Thresholds for viewing isoforms in report table - transcript_table_cov_thresh = 50 + isoform_table_nrows = 5000 out_dir = "output" sample = null diff --git a/nextflow_schema.json b/nextflow_schema.json index 4d0d2c0..da781aa 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -84,10 +84,10 @@ "description": "Extra options for stringtie transcript assembly.", "default": " --conservative " }, - "transcript_table_cov_thresh": { + "isoform_table_nrows": { "type": "integer", - "description": "Minimum coverage for a transcript to appear in the report table", - "default": 50 + "description": "Maximum rows to dispay in the isoform report table", + "default": 5000 }, "denovo": { "type": "boolean",