Switch isoform table thresholding to just displaying to n rows based on coverage

This commit is contained in:
Neil Horner 2022-08-08 13:48:15 +00:00
parent fad574dacf
commit d477a747e6
4 changed files with 29 additions and 44 deletions

View File

@ -6,7 +6,7 @@ from collections import Counter, defaultdict, OrderedDict
import math import math
from pathlib import Path 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 import simple as scomponents
from aplanat.components.fastcat import read_length_plot, read_quality_plot from aplanat.components.fastcat import read_length_plot, read_quality_plot
from aplanat.report import WFReport from aplanat.report import WFReport
@ -594,21 +594,25 @@ def cluster_quality(cluster_qc_dir, report, sample_ids):
section.plot(cover_panel) section.plot(cover_panel)
def transcript_table(report, df_tmaps, covr_threshold): def transcript_table(report, df_tmaps, max_rows):
"""Create searchable table of transcripts.""" """Create searchable table of transcripts.
:param df_tmaps: pd.DataFrame of all gffcomapre `.tmap` files from all
samples
"""
section = report.add_section() 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 # all in single table and sample_id column? Currently it's the latter
# drop some columns for the big table and do some filtering # drop some columns for the big table and do some filtering
section.markdown(''' section.markdown('''
### Isoforms table ### Isoforms table
Low coverage transcripts are removed to speed up the table viewing. <br> Table interactivity can be slow if too many isoforms are loaded. <br>
Coverage threshold can be set with the parameter The number of isoform rows to load in this table can be set with
`transcript_table_cov_thresh`. `isoform_table_nrows`. It is currently set to `{}`
''') '''.format(max_rows))
df = df_tmaps.drop( df = df_tmaps.drop(
columns=[ columns=[
@ -619,31 +623,6 @@ def transcript_table(report, df_tmaps, covr_threshold):
section.markdown("No transcripts found") section.markdown("No transcripts found")
return 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 # Make a column of number of isoforms in parent gene
gb = df.groupby(['ref_gene_id', 'sample_id']).count() gb = df.groupby(['ref_gene_id', 'sample_id']).count()
# gb = gb.set_index(['ref_gene_id', 'sample_id']) # 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 # Uncalssified transcritps should not be lumped togetehr
df.loc[df.class_code == 'u', 'parent gene iso num'] = None 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) section.table(df, index=False)
@ -664,7 +649,7 @@ def transcriptome_summary(report, gffs, sample_ids, denovo=False):
Plot transcriptome summaries. Plot transcriptome summaries.
Some of this data is available via gffcompare output, but the de novo 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 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 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( bar_isos = hist.histogram(
[isoforms_per_gene], colors=[Colors.cerulean], [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.xaxis.axis_label = "Num. isoforms"
bar_isos.yaxis.axis_label = "Num. genes" bar_isos.yaxis.axis_label = "Num. genes"
@ -864,8 +849,8 @@ def main():
"--sample_ids", required=True, nargs='+', "--sample_ids", required=True, nargs='+',
help="List of sample ids") help="List of sample ids")
parser.add_argument( parser.add_argument(
"--transcript_table_cov_thresh", required=False, type=int, default=50, "--isoform_table_nrows", required=False, type=int, default=5000,
help="Isoforms without this support will be excluded from the table") help="Maximum rows to display in isoforms table")
parser.add_argument( parser.add_argument(
"--cluster_qc_dirs", required=False, type=str, default=None, nargs='*', "--cluster_qc_dirs", required=False, type=str, default=None, nargs='*',
help="Directory with various cluster quality csvs") help="Directory with various cluster quality csvs")
@ -912,7 +897,7 @@ def main():
report.write(args.report) report.write(args.report)
if df_tmaps is not None: 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: if args.cluster_qc_dirs is not None:
cluster_quality(args.cluster_qc_dirs, report, sample_ids) cluster_quality(args.cluster_qc_dirs, report, sample_ids)

View File

@ -309,7 +309,7 @@ process makeReport {
--summaries $seq_summaries \ --summaries $seq_summaries \
--gffcompare_dir $gffcmp_dir \ --gffcompare_dir $gffcmp_dir \
--gff_annotation $gff_annotation \ --gff_annotation $gff_annotation \
--transcript_table_cov_thresh $params.transcript_table_cov_thresh \ --isoform_table_nrows $params.isoform_table_nrows \
$OPT_JAFFAL_CSV \ $OPT_JAFFAL_CSV \
$OPT_DENOVO $OPT_DENOVO
""" """

View File

@ -17,7 +17,7 @@ params {
ref_annotation = null ref_annotation = null
threads = 4 threads = 4
// Thresholds for viewing isoforms in report table // Thresholds for viewing isoforms in report table
transcript_table_cov_thresh = 50 isoform_table_nrows = 5000
out_dir = "output" out_dir = "output"
sample = null sample = null

View File

@ -84,10 +84,10 @@
"description": "Extra options for stringtie transcript assembly.", "description": "Extra options for stringtie transcript assembly.",
"default": " --conservative " "default": " --conservative "
}, },
"transcript_table_cov_thresh": { "isoform_table_nrows": {
"type": "integer", "type": "integer",
"description": "Minimum coverage for a transcript to appear in the report table", "description": "Maximum rows to dispay in the isoform report table",
"default": 50 "default": 5000
}, },
"denovo": { "denovo": {
"type": "boolean", "type": "boolean",