diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2e6153c..315a68c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,14 @@ 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.1.0/),
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+## [Unreleased]
+### Updated
+- Workflow report updated to use `ezcharts`.
+### Fixed
+- Exons per isoforms histogram reporting incorrect numbers.
+### Removed
+- per-class gffcompare tracking files as there exists a combine tracking file.
+
## [v1.4.0]
## Added
- `--igv` parameter (default: false) for outputting IGV config allowing visualisation of read alignments in the EPI2ME App.
diff --git a/README.md b/README.md
index c31271a..bd8af7b 100644
--- a/README.md
+++ b/README.md
@@ -210,13 +210,13 @@ Output files may be aggregated including information for all samples or provided
| Title | File path | Description | Per sample or aggregated |
|-------|-----------|-------------|--------------------------|
| workflow report | wf-transcriptomes-report.html | a HTML report document detailing the primary findings of the workflow | aggregated |
-| Per file read stats | fastq_ingress_results/reads/fastcat_stats/per-file-stats.tsv | A TSV with per file read stats, including all samples. | aggregated |
-| Read stats | fastq_ingress_results/reads/fastcat_stats/per-read-stats.tsv | A TSV with per read stats, including all samples. | aggregated |
-| Run ID's | fastq_ingress_results/reads/fastcat_stats/run_ids | List of run IDs present in reads. | aggregated |
-| Meta map json | fastq_ingress_results/reads/metamap.json | Metadata used in workflow presented in a JSON. | aggregated |
-| Concatenated sequence data | fastq_ingress_results/reads/{{ alias }}.fastq.gz | Per sample reads concatenated in to one FASTQ file. | per-sample |
-| Assembled transcriptome | {{ alias }}_transcriptome.fas | Per sample assembled transcriptome. | per-sample |
-| Annotated assembled transcriptome | {{ alias }}_merged_transcriptome.fas | Per sample annotated assembled transcriptome. | per-sample |
+| Per file read stats | fastq_ingress_results/{{ alias }}//reads/fastcat_stats/per-file-stats.tsv | A TSV with per file read stats, including all samples. | aggregated |
+| Read stats | fastq_ingress_results/{{ alias }}//reads/fastcat_stats/per-read-stats.tsv | A TSV with per read stats, including all samples. | aggregated |
+| Run ID's | fastq_ingress_results/{{ alias }}//reads/fastcat_stats/run_ids | List of run IDs present in reads. | aggregated |
+| Meta map json | fastq_ingress_results/{{ alias }}//reads/metamap.json | Metadata used in workflow presented in a JSON. | aggregated |
+| Concatenated sequence data | fastq_ingress_results/{{ alias }}//reads/{{ alias }}.fastq.gz | Per sample reads concatenated in to one FASTQ file. | per-sample |
+| Assembled transcriptome | {{ alias }}_transcriptome.fas | Per sample assembled transcriptome. Not output if a reference annotation was supplied | per-sample |
+| Annotated assembled transcriptome | {{ alias }}_merged_transcriptome.fas | Per sample annotated assembled transcriptome. Only output if a reference annotation was supplied | per-sample |
| Alignment summary statistics | {{ alias }}_read_aln_stats.tsv | Per sample alignment summary statistics. | per-sample |
| GFF compare results. | {{ alias }}_gffcompare | All GFF compare output files. | per-sample |
| Differential gene expression results | de_analysis/results_dge.tsv | This is a gene-level result file that describes genes and their probability of showing differential expression between experimental conditions. | aggregated |
diff --git a/bin/workflow_glue/compute_cluster_quality.py b/bin/workflow_glue/compute_cluster_quality.py
deleted file mode 100755
index 9e768de..0000000
--- a/bin/workflow_glue/compute_cluster_quality.py
+++ /dev/null
@@ -1,649 +0,0 @@
-#!/usr/bin/env python
-"""Generate cluster quality data."""
-
-# Adapted form script by Kristoffer Sahlin for
-# isONclust: https://github.com/ksahlin/isONclust
-
-from collections import defaultdict
-import math
-from pathlib import Path
-
-import matplotlib
-from matplotlib import pyplot as plt
-from matplotlib.backends.backend_pdf import PdfPages
-import pandas as pd
-import pysam
-from sklearn.metrics.cluster import (
- adjusted_rand_score, completeness_score,
- homogeneity_score, v_measure_score)
-
-from .util import wf_parser, get_named_logger # noqa: ABS101
-
-matplotlib.use('Agg')
-
-logger = get_named_logger("clustqual")
-
-
-def argparser():
- """Argument parser for entrypoint."""
- parser = wf_parser("compute_cluster_quality")
- parser.add_argument(
- '--clusters',
- type=str,
- help='Inferred clusters (tsv file)')
- parser.add_argument(
- '--classes',
- type=str,
- help='A sorted and indexed bam file.')
- parser.add_argument(
- '--ctsv',
- default=None,
- type=str,
- help='Write true classes in this TSV file.')
- parser.add_argument(
- '--simulated',
- action="store_true",
- help='Simulated data, we can simply read correct classes '
- 'from the ref field.')
- parser.add_argument(
- '--ont',
- action="store_true",
- help='ONT data, parsing accessions differently.')
- parser.add_argument(
- '--modified_ont',
- action="store_true",
- help='ONT data preprocessed accessions, parsing '
- 'accessions differently.')
- parser.add_argument('--outfile', type=str, help='Output file with results')
- parser.add_argument(
- '--report',
- type=str,
- help='Output PDF file with report')
- parser.add_argument('--sizes', type=str, help='Cluster sizes')
- parser.add_argument(
- '--raw_data_out',
- type=str,
- help='dir to save raw data for plotting')
-
- return parser
-
-
-def parse_inferred_clusters_tsv(tsv_file, args):
- """parse_inferred_clusters_tsv."""
- infile = open(tsv_file, "r")
- infile.readline()
- clusters = {}
- for line in infile:
- cluster_id, _, read_acc = line.strip().split("\t")
- if args.simulated:
- read_acc = "_".join([item for item in read_acc.split("_")[:-1]])
- elif args.ont:
- read_acc = read_acc.split(" ")[0]
- elif args.modified_ont:
- read_acc = read_acc
- else:
- read_acc = read_acc.split("_strand")[0]
-
- clusters[read_acc] = int(cluster_id)
- return clusters
-
-
-def parse_true_clusters(ref_file):
- """parse_true_clusters."""
- classes = defaultdict(dict)
- ref_id_to_chrom = {}
- alignment_counter = defaultdict(int)
-
- prev_chrom = -1
- curr_class_id = -1
- prev_class_stop = -1
- prev_read_id = ""
- unique_reads = set()
- unclassified = 0
- for read in ref_file.fetch(until_eof=True):
- unique_reads.add(read.query_name)
- if read.is_unmapped:
- unclassified += 1
- continue
- # deal with supplementary alignments!!
- if read.is_secondary or read.is_supplementary:
- continue
- assert prev_read_id != read.query_name
-
- chrom = read.reference_name
- if chrom != prev_chrom:
- curr_class_id += 1
- classes[read.query_name] = curr_class_id
- prev_chrom = chrom
- prev_class_stop = read.reference_end
-
- else:
- read_ref_start = read.reference_start
- if read_ref_start > prev_class_stop:
- curr_class_id += 1
- classes[read.query_name] = curr_class_id
- prev_class_stop = read.reference_end
- else:
- classes[read.query_name] = curr_class_id
-
- prev_class_stop = max(read.reference_end, prev_class_stop)
- prev_read_id = read.query_name
-
- # classes[read.query_name] = int(read.reference_id) # chrom
- ref_id_to_chrom[int(read.reference_id)] = chrom
- alignment_counter[int(read.reference_id)] += 1
- # if chrom not in class_ranges:
- # class_ranges[chrom] = {}
-
- # for start, stop in class_ranges[chrom]:
- # if start <= read_ref_start and read_ref_end <= stop:
- # # entirly within
- # elif
-
- # classes[read.query_name] = #read.reference_name.split("|")[0]
- # #"|".join(read.reference_name.split("|")[:2])
- return classes, len(unique_reads), unclassified
-
-
-def parse_true_clusters_simulated(ref_file):
- """parse_true_clusters_simulated."""
- classes = defaultdict(dict)
- for read in ref_file.fetch(until_eof=True):
- # by gene id
- classes[read.query_name] = read.reference_name.split("|")[0]
- # by transcript id
- # classes[read.query_name] = read.reference_name.split("|")[1]
-
- return classes
-
-
-def compute_v_measure(clusters, classes):
- """compute_v_measure."""
- class_list, cluster_list = [], []
- # not_found_id = 1000000
- clustered_but_unaligned = 0
- for read in clusters:
- if read in classes:
- class_list.append(classes[read])
- cluster_list.append(clusters[read])
- else:
- clustered_but_unaligned += 1
-
- # added the unprocessed reads to the measure
- not_clustered = set(classes.keys()) - set(clusters.keys())
- highest_cluster_id = max(clusters.values())
- highest_cluster_id += 1
- for read in not_clustered:
- class_list.append(classes[read])
- cluster_list.append(highest_cluster_id)
- highest_cluster_id += 1
-
- v_score = v_measure_score(class_list, cluster_list)
- compl_score = completeness_score(class_list, cluster_list)
- homog_score = homogeneity_score(class_list, cluster_list)
- ari = adjusted_rand_score(class_list, cluster_list)
- logger.info("Not included in clustering but aligned:", len(not_clustered))
- logger.info(
- "v:",
- v_score,
- "Completeness:",
- compl_score,
- "Homogeneity:",
- homog_score)
- logger.info(
- "Nr reads clustered but unaligned "
- "(i.e., no class and excluded from v-measure): ",
- clustered_but_unaligned)
- return v_score, compl_score, homog_score, clustered_but_unaligned, ari
-
-
-def compute_v_measure_non_singleton_classes(clusters, classes):
- """V measure for non-singleton classes."""
- max_cluster_id = max(clusters.values())
- new_id = max_cluster_id + 1
- classes_dict = {}
- for read_acc, cl_id in classes.items():
- if cl_id not in classes_dict:
- classes_dict[cl_id] = [read_acc]
- else:
- classes_dict[cl_id].append(read_acc)
-
- nontrivial_classes_reads = []
- for cl_id in classes_dict:
- if len(classes_dict[cl_id]) < 5:
- continue
- else:
- for read in classes_dict[cl_id]:
- nontrivial_classes_reads.append(read)
-
- class_list, cluster_list = [], []
- for read in nontrivial_classes_reads:
- if read in clusters:
- class_list.append(classes[read])
- cluster_list.append(clusters[read])
- else:
- class_list.append(classes[read])
- cluster_list.append(new_id)
- new_id += 1
-
- v_score = v_measure_score(class_list, cluster_list)
- compl_score = completeness_score(class_list, cluster_list)
- homog_score = homogeneity_score(class_list, cluster_list)
- nr_filtered_classes = len(
- [1 for cl_id in classes_dict if len(classes_dict[cl_id]) >= 5])
- logger.info(
- "NONTRIvIAL CLASSES: v:",
- v_score,
- "Completeness:",
- compl_score,
- "Homogeneity:",
- homog_score)
- logger.info("NUMBER OF CLASSES (FILTERED):", len(
- [1 for cl_id in classes_dict if len(classes_dict[cl_id]) >= 5]))
- return v_score, compl_score, homog_score, nr_filtered_classes
-
-
-def compute_v_measure_non_singletons(clusters, classes):
- """V measure for non-singletons."""
- cluster_dict = {}
- for read_acc, cl_id in clusters.items():
- if cl_id not in cluster_dict:
- cluster_dict[cl_id] = [read_acc]
- else:
- cluster_dict[cl_id].append(read_acc)
-
- nontrivial_clustered_reads = []
- for cl_id in cluster_dict:
- if len(cluster_dict[cl_id]) <= 1:
- continue
- else:
- for read in cluster_dict[cl_id]:
- nontrivial_clustered_reads.append(read)
-
- class_list, cluster_list = [], []
- # not_found_id = 1000000
- clustered_but_unaligned = 0
- for read in nontrivial_clustered_reads:
- if read in classes:
- class_list.append(classes[read])
- cluster_list.append(clusters[read])
- else:
- logger.debug("Read was clustered but unaligned:", read)
- clustered_but_unaligned += 1
-
- v_score = v_measure_score(class_list, cluster_list)
- compl_score = completeness_score(class_list, cluster_list)
- homog_score = homogeneity_score(class_list, cluster_list)
- logger.info(
- "NONTRIvIAL CLUSTERS: v:",
- v_score,
- "Completeness:",
- compl_score,
- "Homogeneity:",
- homog_score)
- logger.info(
- "NONTRIvIAL CLUSTERS: Nr reads clustered but unaligned "
- "(i.e., no class and excluded from v-veasure): ",
- clustered_but_unaligned)
- return v_score, compl_score, homog_score, clustered_but_unaligned
-
-
-def percentile(n, percent, key=lambda x: x):
- """
- Find the percentile of a list of values.
-
- @parameter n - is a list of values. Note N MUST BE already sorted.
- @parameter percent - a float value from 0.0 to 1.0.
- @parameter key - optional key function to compute value
- from each element of N.
-
- @return - the percentile of the values
- """
- if not n:
- return None
- k = (len(n) - 1) * percent
- f = math.floor(k)
- c = math.ceil(k)
- if f == c:
- return key(n[int(k)])
- d0 = key(n[int(f)]) * (c - k)
- d1 = key(n[int(c)]) * (k - f)
- return d0 + d1
-
-# end of http://code.activestate.com/recipes/511478/ }}}
-
-
-def get_cluster_information(clusters, classes):
- """Get cluster info."""
- # class distribution
- class_dict = {}
- for read_acc, class_id in classes.items():
- if class_id not in class_dict:
- class_dict[class_id] = [read_acc]
- else:
- class_dict[class_id].append(read_acc)
-
- total_nr_classes = len(class_dict)
- class_distribution = sorted([len(cl) for cl in class_dict.values()])
- singleton_classes = set(
- [acc_list[0] for cl_id, acc_list in class_dict.items()
- if len(acc_list) == 1])
- min_class_size = min(class_distribution)
- max_class_size = max(class_distribution)
- mean_class_size = sum(class_distribution) / float(len(class_distribution))
- median_class_size = class_distribution[int(len(class_distribution) / 2)] \
- if len(class_distribution) % 2 == 1 else (
- class_distribution[int(len(class_distribution) / 2)] +
- class_distribution[int(len(class_distribution) / 2) - 1]) / 2.0
-
- upper_75_class_size = percentile(class_distribution, 0.75)
- median_class_size = percentile(class_distribution, 0.5)
-
- tot_size = sum(class_distribution)
- e_class_size = sum(
- [c_s**2 for c_s in class_distribution]) / float(tot_size)
- tot_iterated_size = 0
- for c_s in class_distribution[::-1]:
- tot_iterated_size += c_s
- if tot_iterated_size >= tot_size / 2.0:
- n50_class_size = c_s
- break
-
- # cluster distribution
- cluster_dict = {}
- for read_acc, cl_id in clusters.items():
- if cl_id not in cluster_dict:
- cluster_dict[cl_id] = [read_acc]
- else:
- cluster_dict[cl_id].append(read_acc)
-
- cluster_distribution = sorted([len(cl) for cl in cluster_dict.values()])
-
- # in case unclustered reads are missing from output (as for isoseq3)
- omitted_from_output_singletons = set(classes.keys()) - set(clusters.keys())
- cluster_distribution = [
- 1 for i in range(
- len(omitted_from_output_singletons))] + cluster_distribution
- total_nr_clusters = len(cluster_distribution)
-
- singleton_clusters = set(
- [acc_list[0] for cl_id, acc_list in cluster_dict.items()
- if len(acc_list) == 1])
- min_cluster_size = min(cluster_distribution)
- max_cluster_size = max(cluster_distribution)
- mean_cluster_size = sum(cluster_distribution) / \
- float(len(cluster_distribution))
-
- upper_75_cluster_size = percentile(cluster_distribution, 0.75)
- median_cluster_size = percentile(cluster_distribution, 0.5)
-
- tot_size = sum(cluster_distribution)
- e_cluster_size = sum(
- [c_s**2 for c_s in cluster_distribution]) / float(tot_size)
- tot_iterated_size = 0
- for c_s in cluster_distribution[::-1]:
- tot_iterated_size += c_s
- if tot_iterated_size >= tot_size / 2.0:
- n50_cluster_size = c_s
- break
-
- unaligned_but_nontrivially_clustered = set(
- clusters.keys()) - singleton_clusters - set(classes.keys())
-
- # not_considered = set([read for read in classes if read not in clusters ])
-
- not_clustered_classes = defaultdict(int)
- clustered_classes = defaultdict(int)
- reads_not_clustered = defaultdict(list)
- for read in clusters:
- if read in classes:
- class_id = classes[read]
- else:
- class_id = "unaligned"
-
- if read in singleton_clusters:
- not_clustered_classes[class_id] += 1
- reads_not_clustered[class_id].append(read)
- else:
- clustered_classes[class_id] += 1
-
- logger.info("UNCLUSTERED:", "Tot classes:", len(not_clustered_classes))
- logger.info("CLUSTERED:", "Tot classes:", len(clustered_classes))
- logger.info("MIXED:", "Tot classes containing both:", len(
- set(clustered_classes.keys()) & set(not_clustered_classes.keys())))
- logger.info("Total number of classes (unique gene ID):", total_nr_classes)
- return (
- total_nr_classes - len(singleton_classes),
- len(singleton_classes),
- min_class_size,
- max_class_size,
- mean_class_size,
- median_class_size,
- total_nr_clusters,
- len(singleton_clusters) + len(omitted_from_output_singletons),
- min_cluster_size,
- max_cluster_size,
- mean_cluster_size,
- median_cluster_size,
- len(unaligned_but_nontrivially_clustered),
- upper_75_class_size,
- upper_75_cluster_size,
- e_class_size,
- n50_class_size,
- e_cluster_size,
- n50_cluster_size)
-
-
-def main(args):
- """Entry point."""
- clusters = parse_inferred_clusters_tsv(args.clusters, args)
- if not clusters:
- outfile = open(args.outfile, "w")
- outfile.write("No clusters created\n")
- outfile.close()
- return
-
- if args.simulated:
- ref_file = pysam.AlignmentFile(args.classes, "r", check_sq=False)
- classes = parse_true_clusters_simulated(ref_file)
- # by simulation we know classes of all reads, they are therefore the
- # same number.
- tot_nr_reads = len(classes)
- else:
- ref_file = pysam.AlignmentFile(args.classes, "rb", check_sq=False)
- classes, tot_nr_reads, unclassified = parse_true_clusters(ref_file)
-
- v_score, compl_score, homog_score, clustered_but_unaligned, ari = \
- compute_v_measure(clusters, classes)
-
- (
- nr_non_singleton_classes, singleton_classes, min_class_size,
- max_class_size, mean_class_size, median_class_size, total_nr_clusters,
- singleton_clusters, min_cluster_size, max_cluster_size,
- mean_cluster_size, median_cluster_size,
- unaligned_but_nontrivially_clustered,
- upper_75_class_size, upper_75_cluster_size, e_class_size,
- n50_class_size, e_cluster_size,
- n50_cluster_size
- ) = get_cluster_information(clusters, classes)
-
- outfile = open(args.outfile, "w")
-
- outfile.write("CLASSES\n")
-
- # reads, unaligned, classes, singleton, min,max, mean,median
-
- outfile.write(
- "{0},{1},{2},{3},{4},{5},{6},{7}\n".format(
- "tot_nr_reads",
- "unclassified",
- "nr_non_singleton_classes",
- "singleton_classes",
- "upper_75_class_size",
- "median_class_size",
- "e_class_size",
- "n50_class_size"))
- outfile.write(
- "{0},{1},{2},{3},{4},{5},{6},{7}\n".format(
- tot_nr_reads,
- unclassified,
- nr_non_singleton_classes,
- singleton_classes,
- upper_75_class_size,
- median_class_size,
- e_class_size,
- n50_class_size))
-
- # reads_nontrivially_clustered_(%), Singletons_(%),
- # reads_nontrivially_clustered_but_unaligned, v, c,h ,v_nt, c_nt,h_nt,
- # non_singleton_clusters, min, max, median, mean
-
- reads_nontrivially_clustered_percent = round(
- 100 * (float(tot_nr_reads - singleton_clusters) / tot_nr_reads), 1)
- # round(1.0 - reads_nontrivially_clustered_percent, 2)
- reads_nontrivially_clustered_but_unaligned = \
- unaligned_but_nontrivially_clustered
- v, c, h = round(v_score, 3), round(compl_score, 3), round(homog_score, 3)
-
- non_singleton_clusters = total_nr_clusters - singleton_clusters
-
- logger.info(
- "NONTRIVIAL CLUSTERS: ", (total_nr_clusters - singleton_clusters))
-
- outfile.write("CLUSTERS\n")
- outfile.write(
- "{0},{1},{2},{3},{4},{5},{6},{7},{8},{9},{10},{11}\n".format(
- "v",
- "c",
- "h",
- "ARI",
- "reads_nontrivially_clustered_percent",
- "reads_nontrivially_clustered_but_unaligned",
- "non_singleton_clusters",
- "singleton_clusters",
- "upper_75_cluster_size",
- "median",
- "e_cluster_size",
- "n50_cluster_size"))
- outfile.write(
- "{0},{1},{2},{3},{4},{5},{6},{7},{8},{9},{10},{11}\n".format(
- v,
- c,
- h,
- ari,
- reads_nontrivially_clustered_percent,
- reads_nontrivially_clustered_but_unaligned,
- non_singleton_clusters,
- singleton_clusters,
- upper_75_cluster_size,
- median_cluster_size,
- e_cluster_size,
- n50_cluster_size))
- outfile.close()
-
- if args.ctsv is not None:
- cfh = open(args.ctsv, "w")
- cfh.write("Read\tCluster\n")
- for r, c in classes.items():
- cfh.write("{}\t{}\n".format(r, c))
- cfh.flush()
- cfh.close()
-
- dfc = pd.DataFrame(
- {
- 'Statistic': [
- 'v-measure',
- 'ARI',
- 'Completeness',
- 'Homogeneity'],
- 'value': [
- v,
- ari,
- c,
- h]}).set_index('Statistic')
- dfn = pd.DataFrame(
- {'Statistic': ['NonSingleton',
- 'Singletons'],
- 'value': [non_singleton_clusters,
- singleton_clusters]}).set_index('Statistic')
- dfs = pd.DataFrame(
- {
- 'Statistic': [
- 'Upper75ClsSize',
- 'Upper75ClassSize',
- 'MedianClsSize',
- 'MedianClassSize'],
- 'value': [
- upper_75_cluster_size,
- upper_75_class_size,
- median_cluster_size,
- median_class_size]}).set_index('Statistic')
- dfs2 = pd.DataFrame(
- {'Statistic': ['N50ClsSize', 'N50ClassSize'],
- 'value': [n50_cluster_size, n50_class_size]
- }).set_index('Statistic')
-
- rdo = Path(args.raw_data_out)
- dfc.to_csv(rdo / 'v_ari_com_hom.csv')
- dfn.to_csv(rdo / 'sing_nonsing.csv')
- dfs.to_csv(rdo / 'class_sizes1.csv')
- dfs2.to_csv(rdo / 'class_sizes2.csv')
-
- pages = PdfPages(args.report)
-
- yd = 7 * 2.5
-
- ax = dfc.plot(kind='bar', fontsize=7, rot=0)
- for p in ax.patches:
- ax.annotate(
- "{:.2f}".format(
- p.get_height()),
- (p.get_x() +
- p.get_width() /
- 2.0,
- p.get_height()),
- ha='center')
- pages.savefig()
- plt.clf()
-
- ax = dfn.plot(kind='bar', fontsize=7, rot=0)
- for p in ax.patches:
- ax.annotate(
- "{:.2f}".format(
- p.get_height()),
- (p.get_x() +
- p.get_width() /
- 2.0,
- p.get_height()),
- ha='center')
- pages.savefig()
- plt.clf()
-
- ax = dfs.plot(kind='bar', fontsize=7, rot=0)
- for p in ax.patches:
- ax.annotate(
- "{:.2f}".format(
- p.get_height()),
- (p.get_x() +
- p.get_width() /
- 2.0,
- p.get_height() +
- yd),
- ha='center')
- pages.savefig()
- plt.clf()
- ax = dfs2.plot(kind='bar', fontsize=7, rot=0)
- for p in ax.patches:
- ax.annotate(
- "{:.2f}".format(
- p.get_height()),
- (p.get_x() +
- p.get_width() /
- 2.0,
- p.get_height() +
- yd),
- ha='center')
- pages.savefig()
- plt.clf()
-
- pages.close()
diff --git a/bin/workflow_glue/configure_igv.py b/bin/workflow_glue/configure_igv.py
index db158a4..db735dd 100755
--- a/bin/workflow_glue/configure_igv.py
+++ b/bin/workflow_glue/configure_igv.py
@@ -1,135 +1,313 @@
"""Create an IGV config file."""
-from itertools import zip_longest
import json
+from pathlib import Path
import sys
from .util import get_named_logger, wf_parser # noqa: ABS101
-def parse_fnames(fofn):
- """Parse list with filenames and return them grouped as ref-, XAM-, or VCF-related.
+# Common variables
+REF_EXTENSIONS = [".fasta", ".fasta.gz", ".fa", ".fa.gz", ".fna", ".fna.gz"]
+DATA_TYPES_LISTS = {
+ "bam": ["bam"],
+ "bam_idx": ["bam.bai"],
+ "cram": ["cram"],
+ "cram_idx": ["cram.crai"],
+ "vcf": ["vcf", "vcf.gz"],
+ "vcf_idx": ["vcf.gz.tbi", "vcf.gz.csi"],
+ "bcf": ["bcf"],
+ "bcf_idx": ["bcf.csi"],
+ "gtf": ["gtf", "gtf.gz"],
+ "gtf_idx": ["gtf.gz.tbi"],
+ "gff": ["gff", "gff.gz", "gff3", "gff3.gz"],
+ "gff_idx": ["gff.gz.tbi", "gff3.gz.tbi"],
+ "bed": ["bed", "bed.gz"],
+ "bed_idx": ["bed.gz.tbi"],
+ "bedmethyl": ["bedmethyl", "bedmethyl.gz"],
+ "bedmethyl_idx": ["bedmethyl.gz.tbi"],
+ "ref": REF_EXTENSIONS,
+}
+DATA_TYPES = {
+ ext: ftype for ftype, extlist in DATA_TYPES_LISTS.items() for ext in extlist
+}
- :param fofn: File with list of file names (one per line)
- :return: dict of reference-related filenames (with keys 'ref', 'fai', and '.gzi' and
- `None` as default values); lists of XAM- and VCF-related filenames
- """
- ref_extensions = [".fasta", ".fasta.gz", ".fa", ".fa.gz", ".fna", ".fna.gz"]
- ref_dict = {}
- xams = []
- xam_indices = []
- vcfs = []
- vcf_indices = []
- with open(fofn, "r") as f:
- for line in f:
- fname = line.strip()
- if any(fname.endswith(ext) for ext in ref_extensions):
- ref_dict["ref"] = fname
- elif fname.endswith(".fai"):
- ref_dict["fai"] = fname
- elif fname.endswith(".gzi"):
- ref_dict["gzi"] = fname
- elif fname.endswith(".bam") or fname.endswith(".cram"):
- xams.append(fname)
- elif fname.endswith(".bai") or fname.endswith(".crai"):
- xam_indices.append(fname)
- elif fname.endswith(".vcf") or fname.endswith(".vcf.gz"):
- vcfs.append(fname)
- elif fname.endswith(".csi") or fname.endswith(".tbi"):
- vcf_indices.append(fname)
- # do some sanity checks
- if "ref" not in ref_dict:
- raise ValueError(
- "No reference file (i.e. file ending in one of "
- f"{ref_extensions} was found)."
- )
- ref = ref_dict["ref"]
- if (gzi := ref_dict.get("gzi")) is not None:
- # since we got a '.gzi' index, make sure that the reference is actually
- # compressed
- if not ref_dict["ref"].endswith(".gz"):
- raise ValueError(
- f"Found GZI reference index '{gzi}', but the reference file "
- f"'{ref}' appears not to be compressed."
+# Data by idx
+DATA_INDEXES_FMT = {
+ fmt: f"{fmt}_idx" for fmt, dtype in DATA_TYPES.items() if "_idx" not in dtype
+}
+
+# Assign each format to its index
+INDEX_PAIRS = {
+ "bam": ("bai",),
+ "cram": ("crai",),
+ "vcf": ("tbi", "csi"),
+ "bcf": ("csi",),
+ "bed": ("tbi",),
+ "bedmethyl": ("tbi",),
+ "gff": ("tbi",),
+ "gtf": ("tbi",),
+}
+
+
+class TrackBuilder:
+ """Class that builds an IGV track."""
+
+ def __init__(self):
+ """Initialize properties for interval track."""
+ # Reference properties
+ self.ref = None
+ self.fai = None
+ self.gzi = None
+ # Samples info
+ self.samples = {}
+ # Track properties
+ self.igv_json = {"reference": {}, "tracks": []}
+ self.track_type = {
+ "bam": "alignment",
+ "cram": "alignment",
+ "bcf": "variant",
+ "vcf": "variant",
+ "bedmethyl": "annotation",
+ "bed": "annotation",
+ "gtf": "annotation",
+ "gff": "annotation",
+ }
+ # Here we save aliases of file formats that IGV.js
+ # wants and that do not match the input file extension.
+ self.igv_fmt_alias = {"gff": "gff3"}
+ # lookup of extra options for each data type
+ self.extra_opts_lookups = {
+ "bam": {},
+ "cram": {},
+ "bcf": {},
+ "vcf": {},
+ "bed": {},
+ "bedmethyl": {},
+ "gtf": {},
+ "gff": {},
+ }
+
+ def add_ref(self, ref=None):
+ """Add reference file, unless already defined."""
+ if self.ref:
+ raise Exception(
+ f"Reference genome has already been set to {self.ref}.\n"
+ "Only one reference FASTA file is expected."
)
- if xam_indices:
- if len(xams) != len(xam_indices):
- raise ValueError("Got different number of XAM and XAM index files.")
- if vcf_indices:
- if len(vcfs) != len(vcf_indices):
- raise ValueError("Got different number of VCF and VCF index files.")
- if xams and vcfs:
- if len(xams) != len(vcfs):
- raise ValueError("Got different number of XAM and VCF files.")
- # if we got XAM or VCF indices, pair them up with their corresponding files (and
- # otherwise with `None`)
- xams_with_indices = zip_longest(xams, xam_indices)
- vcfs_with_indices = zip_longest(vcfs, vcf_indices)
- return ref_dict, xams_with_indices, vcfs_with_indices
+ else:
+ self.ref = ref
+
+ def add_ref_index(self, ref_index=None):
+ """Add reference index if valid."""
+ basename = Path(self.ref).name
+ idx_basename = Path(ref_index).name
+ if idx_basename == f"{basename}.fai":
+ self.fai = ref_index
+ if idx_basename == f"{basename}.gzi" and basename.endswith(".gz"):
+ self.gzi = ref_index
+
+ def parse_fnames(self, fofn):
+ """Parse list with filenames and return them grouped.
+
+ :param fofn: File with list of file names (one per line)
+ """
+ tmp_samples = {}
+ with open(fofn, "r") as f:
+ for line in f:
+ # If the line contains the sample name, prepare the data structure
+ if "," in line:
+ sample, fname = line.strip().split(",")
+ if sample not in tmp_samples:
+ tmp_samples[sample] = SampleBundle(sample=sample)
+ tmp_samples[sample].append(fname)
+ else:
+ # Otherwise, assign everything to NO_SAMPLE
+ # Files will still be displayed, but in no specific order.
+ fname = line.strip()
+ if any(fname.endswith(ext) for ext in REF_EXTENSIONS):
+ self.add_ref(ref=fname)
+ elif fname.endswith(".fai") or fname.endswith(".gzi"):
+ self.add_ref_index(ref_index=fname)
+ else:
+ if "NO_SAMPLE" not in tmp_samples.keys():
+ tmp_samples["NO_SAMPLE"] = SampleBundle(sample="NO_SAMPLE")
+ tmp_samples["NO_SAMPLE"].append(fname)
+ # Re-order samples in dict and add them to the list, leaving
+ # NO_SAMPLE as last
+ sorted_samples = (
+ sorted([sample for sample in tmp_samples.keys() if sample != 'NO_SAMPLE'])
+ )
+ if 'NO_SAMPLE' in tmp_samples.keys():
+ sorted_samples += ['NO_SAMPLE']
+ for sample in sorted_samples:
+ self.samples[sample] = tmp_samples[sample]
+
+ def build_igv_json(self):
+ """Ensure there is a reference genome."""
+ if not self.ref:
+ raise ValueError(
+ "No reference file (i.e. file ending in one of "
+ f"{REF_EXTENSIONS} was found)."
+ )
+ # Evaluate that a bgzipped reference has the appropriate index.
+ if self.ref.endswith(".gz") and not self.gzi:
+ raise ValueError(f"GZI reference index for {self.ref} not found.")
+
+ # Create the base track if there is a reference genome.
+ self.igv_json["reference"] = {
+ "id": "ref",
+ "name": "ref",
+ "wholeGenomeView": False,
+ "fastaURL": self.ref,
+ }
+ if self.fai:
+ self.igv_json["reference"]["indexURL"] = self.fai
+ if self.gzi:
+ self.igv_json["reference"]["compressedIndexURL"] = self.gzi
+
+ # Add samples data now
+ for sample, bundle in self.samples.items():
+ bundle.process_data()
+ # Add the bundled data to the tracks
+ for fname, index, file_fmt in bundle.data_bundles:
+ self.add_track(
+ fname,
+ file_fmt,
+ sample_name=sample if sample != "NO_SAMPLE" else None,
+ index=index,
+ extra_opts=self.extra_opts_lookups[file_fmt],
+ )
+
+ def add_track(self, infile, file_fmt, sample_name=None, index=None, extra_opts={}):
+ """Add a track to an IGV json.
+
+ This function takes an input file, an optional index file, its
+ file format and additional extra options for the track.
+
+ :param infile: input file to create a track for
+ :param file_fmt: input file track type
+ :param sample_name: Name of the sample to display in the track name
+ :param index: index for the input file
+ :param extra_opts: dict of extra options for the track
+ :return: dict with track options
+ """
+ # Define track name depending on whether the sample ID is provided
+ track_name = Path(infile).name
+ if sample_name:
+ track_name = f"{sample_name}: {Path(infile).name}"
+ track_dict = {
+ "name": track_name,
+ "type": self.track_type[file_fmt],
+ "format": self.igv_fmt_alias.get(file_fmt, file_fmt),
+ "url": infile,
+ }
+ # add the index, if present
+ if index:
+ track_dict["indexURL"] = index
+ track_dict.update(extra_opts)
+ self.igv_json["tracks"] += [track_dict]
+
+ def add_locus(self, locus):
+ """Add target locus to the json."""
+ self.igv_json["locus"] = locus
+
+ def add_extra_opts(
+ self,
+ extra_alignment_opts=None,
+ extra_variant_opts=None,
+ extra_interval_opts=None,
+ ):
+ """Import extra options from json files."""
+ if extra_alignment_opts is not None:
+ with open(extra_alignment_opts, "r") as f:
+ extra_alignment_opts_json = json.load(f)
+ for ftype in ["bam", "cram"]:
+ self.extra_opts_lookups[ftype] = extra_alignment_opts_json
+ if extra_variant_opts is not None:
+ with open(extra_variant_opts, "r") as f:
+ extra_variant_opts_json = json.load(f)
+ for ftype in ["vcf", "bcf"]:
+ self.extra_opts_lookups[ftype] = extra_variant_opts_json
+ if extra_interval_opts is not None:
+ with open(extra_interval_opts, "r") as f:
+ extra_interval_opts_json = json.load(f)
+ for ftype in ["bed", "bedmethyl", "gff", "gtf"]:
+ self.extra_opts_lookups[ftype] = extra_interval_opts_json
-def get_reference_options(ref, fai=None, gzi=None):
- """Create dict with IGV reference options.
+class SampleBundle:
+ """Sample data class.
- :param ref: reference file name
- :param fai: name reference `.fai` index file
- :param gzi: name of `.gzi` index file for a compressed reference
- :return: dict with reference options
+ This class stores the data for multiple tracks for a
+ single sample, then is used to generate a collection of
+ IGV.js tracks.
"""
- # initialise the options dict and add the index attributes later
- ref_opts = {
- "id": "ref",
- "name": "ref",
- "wholeGenomeView": False,
- "fastaURL": ref,
- }
- if fai is not None:
- ref_opts["indexURL"] = fai
- if gzi is not None:
- ref_opts["compressedIndexURL"] = gzi
- return ref_opts
+ def __init__(self, sample):
+ """Initialize properties for a sample."""
+ self.sample = sample
+ self.infiles = []
+ self.data_bundles = []
-def get_alignment_track(xam, xai=None, extra_opts=None):
- """Create dict with options for IGV alignment track.
+ def append(self, fname):
+ """Add a new raw file to the bundle."""
+ self.infiles.append(fname)
- :param xam: name of XAM file to be displayed
- :param xai: name of XAM index file
- :param extra_opts: dict of extra options for the alignment track
- :return: dict with alignment track options
- """
- alignment_track_dict = {
- "name": xam,
- "type": "alignment",
- "format": xam.split(".")[-1],
- "url": xam,
- }
- # add the XAM index if present
- if xai is not None:
- alignment_track_dict["indexURL"] = xai
- alignment_track_dict.update(extra_opts or {})
- return alignment_track_dict
+ def process_data(self):
+ """Process input files."""
+ fbasenames = [Path(fname).name for fname in self.infiles]
+ ftypes = [self.classify_files(bname) for bname in fbasenames]
+ self.data_bundles = self.pair_file_with_index(self.infiles, fbasenames, ftypes)
+ @staticmethod
+ def classify_files(fname):
+ """Classify inputs."""
+ for extension, ftype in DATA_TYPES.items():
+ if fname.endswith(f".{extension}"):
+ return ftype
-def get_variant_track(vcf, index=None, extra_opts=None):
- """Create dict with options for IGV variant track.
+ @staticmethod
+ def pair_file_with_index(infiles, fbasenames, ftypes):
+ """Clump files with their indexes."""
+ # Collect data by group type
+ groups = {ftype: {"basenames": [], "paths": []} for ftype in set(ftypes)}
+ # Group each file by its type and base name
+ for ftype, fbasename, fname in zip(ftypes, fbasenames, infiles):
+ groups[ftype]["basenames"] += [fbasename]
+ groups[ftype]["paths"] += [fname]
- :param vcf: name of VCF file to be displayed
- :param index: name of VCF index file (ending in `.csi` or `.tbi`)
- :param extra_opts: dict of extra options for the variant track
- :return: dict with variant track options
- """
- variant_track_dict = {
- "name": vcf,
- "type": "variant",
- "format": "vcf",
- "url": vcf,
- }
- # add the VCF index if we got an index extension
- if index is not None:
- variant_track_dict["indexURL"] = index
- variant_track_dict.update(extra_opts or {})
- return variant_track_dict
+ # Output bundles
+ outputs = []
+ # Start matching the variant files
+ for ftype, itype in DATA_INDEXES_FMT.items():
+ # Ignore file formats that are not present in the bundle.
+ if ftype not in groups:
+ continue
+ # Make pairs of files.
+ for fbasename, fpath in zip(
+ groups[ftype]["basenames"], groups[ftype]["paths"]
+ ):
+ # Construct potential index file names based on basename of input files
+ idx_basenames = set(
+ [f"{fbasename}.{idx}" for idx in INDEX_PAIRS[ftype]]
+ )
+ # Find which indexes are available
+ if itype in groups.keys():
+ idx_basenames = list(
+ idx_basenames.intersection(set(groups[itype]["basenames"]))
+ )
+ # Get the first index (if there are more than one,
+ # it doesn't matter)
+ bname = idx_basenames[0]
+ idx_fn = groups[itype]["paths"][
+ groups[itype]["basenames"].index(bname)
+ ]
+ outputs.append([fpath, idx_fn, ftype])
+ # Otherwise, return only the simple file.
+ else:
+ outputs.append([fpath, None, ftype])
+ return outputs
def main(args):
@@ -137,44 +315,26 @@ def main(args):
logger = get_named_logger("configIGV")
# parse the FOFN
- ref_dict, xams_with_indices, vcfs_with_indices = parse_fnames(args.fofn)
+ igv_builder = TrackBuilder()
+
+ # Add the additional track configurations
+ igv_builder.add_extra_opts(
+ extra_alignment_opts=args.extra_alignment_opts,
+ extra_variant_opts=args.extra_variant_opts,
+ extra_interval_opts=args.extra_interval_opts
+ )
+
+ # Import files
+ igv_builder.parse_fnames(args.fofn)
# initialise the IGV options dict with the reference options
- json_dict = {"reference": get_reference_options(**ref_dict)}
-
- # if we got JSON files with extra options for the alignment / variant tracks, read
- # them
- extra_alignment_opts = {}
- if args.extra_alignment_opts is not None:
- with open(args.extra_alignment_opts, "r") as f:
- extra_alignment_opts = json.load(f)
- extra_variant_opts = {}
- if args.extra_variant_opts is not None:
- with open(args.extra_variant_opts, "r") as f:
- extra_variant_opts = json.load(f)
-
- # now add the alignment and variant tracks
- json_dict["tracks"] = []
- # we use `zip_longest` to make sure that variant and alignment tracks from the same
- # sample are added after each other
- for (vcf, vcf_index), (xam, xam_index) in zip_longest(
- vcfs_with_indices, xams_with_indices, fillvalue=(None, None)
- ):
- if vcf is not None:
- # add a variant track for the VCF
- json_dict["tracks"].append(
- get_variant_track(vcf, vcf_index, extra_variant_opts)
- )
- if xam is not None:
- # add an alignment track for the XAM
- json_dict["tracks"].append(
- get_alignment_track(xam, xam_index, extra_alignment_opts)
- )
+ igv_builder.build_igv_json()
+ # Add locus information
if args.locus is not None:
- json_dict["locus"] = args.locus
+ igv_builder.add_locus(args.locus)
- json.dump(json_dict, sys.stdout, indent=4)
+ json.dump(igv_builder.igv_json, sys.stdout, indent=4)
logger.info("Printed IGV config JSON to STDOUT.")
@@ -202,4 +362,8 @@ def argparser():
"--extra-variant-opts",
help="JSON file with extra options for variant tracks",
)
+ parser.add_argument(
+ "--extra_interval_opts",
+ help="JSON file with extra options for interval tracks",
+ )
return parser
diff --git a/bin/workflow_glue/de_plots.py b/bin/workflow_glue/de_plots.py
index 65c7720..458e7c4 100755
--- a/bin/workflow_glue/de_plots.py
+++ b/bin/workflow_glue/de_plots.py
@@ -1,12 +1,13 @@
#!/usr/bin/env python
"""Create de report section."""
-from glob import glob
import os
-from aplanat import hist, points
-from aplanat.bars import boxplot_series
-from aplanat.util import Colors
+from dominate.tags import h5, p
+from dominate.util import raw
+from ezcharts import scatterplot
+from ezcharts.components.ezchart import EZChart
+from ezcharts.layout.snippets import DataTable
import numpy as np
import pandas as pd
@@ -51,174 +52,55 @@ def create_summary_table(df):
avg_acc, avg_mapq])
-def dtu_table(gene_id, dtu_file, alignment_stats, condition_sheet):
- """Create DTU table and plot."""
- dtu_results = pd.read_csv(dtu_file, sep='\t')
- f = open(condition_sheet)
- df = pd.read_csv(f, sep='\t')
- treated_df = df.loc[df['condition'] == "treated"]
- untreated_df = df.loc[df['condition'] == "untreated"]
- treated_samples = treated_df['sample'].tolist()
- control_samples = untreated_df['sample'].tolist()
- # parmaterise these
- control_name = "condition1"
- treated_name = "condition2"
- table = dtu_results.loc[(dtu_results["geneID"] == gene_id)]
- msg = "Gene ID \"{}\" does not exist in the dataset, please select another"
- assert not table.empty, msg.format(gene_id)
- alignment_stats[alignment_stats["Ref"].isin([gene_id])]
- alignment_stats["Transcript"] = alignment_stats["Ref"].apply(
- lambda x: x.split(".")[0])
- gene_alignments = alignment_stats[alignment_stats["Transcript"].isin(
- table["txID"])]
- gene_alignments = gene_alignments.loc[(
- gene_alignments["Type"] == "Primary")]
- gene_alignments["condition"] = gene_alignments.apply(
- lambda x: control_name if x["fname"] in
- control_samples else treated_name, axis=1)
- groups = gene_alignments.groupby(
- ["Transcript", "fname", "condition"]).agg(
- **{"Transcript_count": ("Read", "size")})
- df = None
- for gene_id, group in gene_alignments.groupby("Transcript"):
- temp = group.groupby("fname").agg(**{
- gene_id: ("Read", "size")
- }).transpose()
- if df is None:
- df = temp
- else:
- df = pd.concat([df, temp])
- df = df.reset_index().rename(columns={"index": "Transcript_ID"})
- table = table.rename(columns={
- "txID": "Transcript_ID",
- "gene": "p_gene",
- "transcript": "p_transcript"
- })
- gene_table = pd.merge(df, table)
- gene_table = gene_table.set_index(["geneID", "Transcript_ID"])
- file_names = set(control_samples.keys())
- file_names.update(treated_samples.keys())
- gene_table = gene_table.fillna(0)
- table = groups.reset_index()[
- ["condition", "Transcript", "Transcript_count"]]
- table['transcript, condition'] = \
- table['Transcript'].astype(str) + ', ' + table['condition'].astype(str)
- repeats = groups.groupby(level=['Transcript', 'condition']).size()
- min_rep, max_rep = min(repeats), max(repeats)
- plot = boxplot_series(
- table['transcript, condition'], table['Transcript_count'],
- x_axis_label='Transcript, condition',
- y_axis_label='Transcript count',
- height=200, width=200,
- title="Transcript counts (from {}-{} replicates)".format(
- min_rep, max_rep))
- plot.xaxis.major_label_orientation = 3.1452/2
- if min_rep < 7:
- for renderer in plot.renderers:
- renderer.glyph.line_alpha = 0.2
- try:
- renderer.glyph.fill_alpha = 0.2
- except Exception:
- pass
- plot.circle(
- table['transcript, condition'], table['Transcript_count'],
- fill_color='black', line_color='black')
- return (table, plot)
-
-
-def pool_csvs(folder):
- """Concat seqkit stats."""
- files = glob(folder + "/*.seqkit.stats")
- dfs = [parse_seqkit(f) for f in files]
- return pd.concat(dfs)
-
-
-def abundance_histogram(filtered_counts, gene_counts, section):
- """Create plot for abundance of transcripts across all samples."""
- section.markdown("""
-Histogram showing the abundance of transcript
-counts for genes identified in the analysis.
- """)
- filtered_count_file = filtered_counts
- gene_count_file = gene_counts
- transcripts_per_gene = pd.read_csv(
- filtered_count_file, sep='\t',
- usecols=['gene_id', 'feature_id']).groupby(['gene_id']).agg(['count'])
- transcripts_per_gene.columns = transcripts_per_gene.columns.droplevel()
- gene_ids = pd.read_csv(
- gene_count_file, sep='\t',
- usecols=[0]).index.values.tolist()
- singletons = [
- gene_id for gene_id in gene_ids if gene_id not in transcripts_per_gene.index.values.tolist()] # noqa
- singletons = pd.DataFrame(index=singletons, columns=['count']).fillna(1)
- transcripts_per_gene = pd.concat([singletons, transcripts_per_gene])
- transcript_plot = hist.histogram(
- [transcripts_per_gene['count'].tolist()],
- binwidth=1, colors=[Colors.cerulean])
- transcript_plot.xaxis.axis_label = "Number of isoforms per gene (n)"
- transcript_plot.yaxis.axis_label = "Number of occurences"
- section.markdown("### Transcripts per gene")
- section.plot(transcript_plot)
-
-
-def dexseq_section(dexseq_file, section, id_dic):
+def dexseq_section(dexseq_file, id_dic, pval_thresh):
"""Add gene isoforms table and plot."""
- section.markdown("### Differential Isoform usage")
- dexseq_caption = '''Table showing gene isoforms, ranked by adjusted
+ h5("Differential Isoform usage")
+ p("""Table showing gene isoforms, ranked by adjusted
p-value, from the DEXSeq analysis. Information shown includes the log2 fold
change between experimental conditions, the log-scaled transcript
- abundance and the false discovery corrected p-value (FDR).
+ abundance and the false discovery corrected p-value (FDR - Benjamini-Hochberg) .
This table has not been filtered
- for genes that satisfy statistical or magnitudinal thresholds'''
- section.markdown(dexseq_caption)
+ for genes that satisfy statistical or magnitudinal thresholds""")
+
dexseq_results = pd.read_csv(dexseq_file, sep='\t')
- dexseq_results.index.name = "gene_id:trancript_id"
+ dexseq_results.index.name = "gene_id:transcript_id"
# Replace gene id with more useful gene name where possible
dexseq_results.index = dexseq_results.index.map(
lambda x: str(id_dic.get(x.split(':')[0])) + ':' + str(x.split(':')[1]))
- dexseq_pvals = dexseq_results.sort_values(by='pvalue', ascending=True)
- section.table(dexseq_results.loc[dexseq_pvals.index], index=True)
- section.markdown("""
-The figure below presents the MA plot from the DEXSeq analysis.
-M is the log2 ratio of isoform transcript abundance between conditions.
-A is the log2 transformed mean abundance value.
-Transcripts that satisfy the logFC and FDR corrected p-value
-thresholds defined are shaded as 'Up-' or 'Down-' regulated.""")
- pval_limit = 0.01
- up = dexseq_results.loc[
- (dexseq_results["Log2FC"] > 0) & (
- dexseq_results['pvalue'] < pval_limit)]
- down = dexseq_results.loc[
- (dexseq_results["Log2FC"] <= 0) & (
- dexseq_results['pvalue'] < pval_limit)]
- not_sig = dexseq_results.loc[(dexseq_results["pvalue"] >= pval_limit)]
- dexseq_plot = points.points(
- x_datas=[
- up["Log2MeanExon"],
- down["Log2MeanExon"],
- not_sig["Log2MeanExon"],
- ],
- y_datas=[
- up["Log2FC"],
- down["Log2FC"],
- not_sig["Log2FC"],
- ],
- title="Average copy per million (CPM) vs Log-fold change (LFC)",
- colors=["red", "blue", "black"],
- names=["Up", "Down", "NotSig"]
- )
+ DataTable.from_pandas(
+ dexseq_results.sort_values(by='pvalue', ascending=True), use_index=True)
- dexseq_plot.xaxis.axis_label = "A (log2 transformed mean exon read counts)"
- dexseq_plot.yaxis.axis_label = """
- M (log2 transformed differential abundance)
- """
- dexseq_results_caption = "### Dexseq results"
- section.markdown(dexseq_results_caption)
- section.plot(dexseq_plot)
+ p(
+ """The figure below presents the MA plot from the DEXSeq analysis.
+ M is the log2 ratio of isoform transcript abundance between conditions.
+ A is the log2 transformed mean abundance value.
+ Transcripts that satisfy the logFC and FDR-corrected
+ (False discovery rate - Benjamini-Hochberg) p-value
+ thresholds defined are shaded as 'Up-' or 'Down-' regulated.""")
+
+ dexseq_results['direction'] = 'not_sig'
+
+ dexseq_results.loc[
+ (dexseq_results["Log2FC"] > 0) & (dexseq_results['pvalue'] < pval_thresh),
+ 'direction'] = 'up'
+
+ dexseq_results.loc[
+ (dexseq_results["Log2FC"] <= 0) & (dexseq_results['pvalue'] < pval_thresh),
+ 'direction'] = 'down'
+
+ plot = scatterplot(
+ data=dexseq_results, x='Log2MeanExon', y='Log2FC', hue='direction',
+ palette=['#E32636', '#7E8896', '#0A22DE'],
+ hue_order=['up', 'down', 'not_sig'], marker='circle')
+ plot._fig.xaxis.axis_label = "A (log2 transformed mean exon read counts)"
+ plot._fig.yaxis.axis_label = "M (log2 transformed differential abundance)"
+ plot.legend = dict(orient='horizontal', top=30)
+ plot._fig.title = "Average copy per million (CPM) vs Log-fold change (LFC)"
+ EZChart(plot)
-def dtu_section(dtu_file, section, gt_dic, ge_dic):
+def dtu_section(dtu_file, gt_dic, ge_dic):
"""Plot dtu section."""
dtu_results = pd.read_csv(dtu_file, sep='\t')
dtu_results["gene_name"] = dtu_results["txID"].apply(
@@ -226,198 +108,187 @@ def dtu_section(dtu_file, section, gt_dic, ge_dic):
dtu_results["geneID"] = dtu_results["geneID"].apply(
lambda x: ge_dic.get(x))
dtu_pvals = dtu_results.sort_values(by='gene', ascending=True)
- dtu_caption = '''Table showing gene and transcript identifiers
- and their FDR corrected probabilities
+ raw("""Table showing gene and transcript identifiers
+ and their FDR-corrected (False discovery rate - Benjamini-Hochberg) probabilities
for the genes and their isoforms that have been
identified as showing DTU using the R packages DEXSeq and StageR.
This list has been shortened requiring that both gene and transcript
must satisfy the p-value
- threshold'''
- section.markdown(dtu_caption)
- section.table(dtu_results.loc[dtu_pvals.index])
+ threshold""")
+ DataTable.from_pandas(dtu_results.loc[dtu_pvals.index], use_index=False)
+
+ raw("""View dtu_plots.pdf file to see plots of differential isoform usage""")
-def copy_dge_add_gene_names(dge_file, geid_gname, newfile_name):
- """Add gene name column to DGE TSV with gene ID to gene Name dict."""
+def dge_section(dge_file, ids_dic, pval_thresh):
+ """Create DGE table and MA plot."""
+ h5("Differential gene expression")
dge_results = pd.read_csv(dge_file, sep='\t')
- column_to_move = dge_results.index.map(
- lambda x: geid_gname.get(x))
- dge_results.insert(0, "gene_name", column_to_move)
- dge_results.to_csv(newfile_name, index=True, index_label="gene_id", sep="\t")
-
-
-def copy_filtered_add_gene_names(unfiltered_file, geid_gname, newfile_name):
- """Use gene id to gene name dict to add name column to counts TSV."""
- unfiltered = pd.read_csv(unfiltered_file, sep='\t')
- unfiltered.insert(1, "gene_name", unfiltered.gene_id.map(
- lambda x: geid_gname.get(x)))
- unfiltered.to_csv(newfile_name, index=False, sep='\t')
-
-
-def copy_tpm_add_gene_names(tpm_file, txid_gname, newfile_name):
- """Use transcript id to gene name dict to add name column to TPM TSV."""
- tpm = pd.read_csv(tpm_file, sep='\t')
- tpm.insert(1, "gene_name", tpm.Reference.map(
- lambda x: txid_gname.get(x)))
- tpm.to_csv(newfile_name, index=False, sep='\t')
-
-
-def dge_section(dge_file, section, ids_dic):
- """Create DGE table and plot."""
- section.markdown('### Differential gene expression')
- dge_results = pd.read_csv(dge_file, sep='\t')
- dge_pvals = dge_results.sort_values(by='FDR', ascending=True)
dge_results[['logFC', 'logCPM', 'F']] = dge_results[
['logFC', 'logCPM', 'F']].round(2)
- dge_caption = """
-Table showing the genes from the edgeR analysis.
-Information shown includes the log2 fold change between
-experimental conditions, the log-scaled counts per million measure of abundance
-and the false discovery corrected p-value (FDR). This table has not been
-filtered for genes that satisfy statistical or magnitudinal thresholds"""
- section.markdown(dge_caption)
+
+ p("""Table showing the genes from the edgeR analysis.
+ Information shown includes the log2 fold change between
+ experimental conditions, the log-scaled counts per million measure of abundance
+ and the FDR-corrected p-value (False discovery rate - Benjamini-Hochberg).
+ This table has not been
+ filtered for genes that satisfy statistical or magnitudinal thresholds""")
+
dge_results.index = dge_results.index.map(lambda x: ids_dic.get(x))
- dge_pvals.index = dge_pvals.index.map(lambda x: ids_dic.get(x))
- section.table(dge_results.loc[dge_pvals.index], index=True)
- dge = pd.read_csv(dge_file, sep="\t")
- section.markdown("""
-This plot visualises differences in measurements between the
-two experimental conditions. M is the log2 ratio of gene expression
-calculated between the conditions.
-A is a log2 transformed mean expression value.
-The figure below presents the MA figure from this edgeR analysis.
-Genes that satisfy the logFC and FDR corrected p-value thresholds
-defined are shaded as 'Up-' or 'Down-' regulated.
+ dge_results = dge_results.sort_values('FDR', ascending=True)
+ dge_results.index.name = 'Transcript'
+ DataTable.from_pandas(dge_results, use_index=True)
+
+ h5("Results of the edgeR Analysis.")
+
+ p("""This plot visualises differences in measurements between the
+ two experimental conditions. M is the log2 ratio of gene expression
+ calculated between the conditions.
+ A is a log2 transformed mean expression value.
+ The figure below presents the MA figure from this edgeR analysis.
+ Genes that satisfy the logFC and FDR-corrected
+ (False discovery rate - Benjamini-Hochberg) p-value thresholds
+ defined are shaded as 'Up-' or 'Down-' regulated.
""")
- pval_limit = 0.01
- up = dge.loc[(dge["logFC"] > 0) & (dge['PValue'] < pval_limit)]
- down = dge.loc[(dge["logFC"] <= 0) & (dge['PValue'] < pval_limit)]
- not_sig = dge.loc[(dge["PValue"] >= pval_limit)]
- logcpm_vs_logfc = points.points(
- x_datas=[
- up["logCPM"],
- down["logCPM"],
- not_sig["logCPM"],
- ],
- y_datas=[
- up["logFC"],
- down["logFC"],
- not_sig["logFC"],
- ],
- title="Average copy per million (CPM) vs Log-fold change (LFC)",
- colors=["red", "blue", "black"],
- names=["Up", "Down", "NotSig"]
- )
- logcpm_vs_logfc.xaxis.axis_label = "Average log CPM"
- logcpm_vs_logfc.yaxis.axis_label = "Log-fold change"
- logcpm_caption = """### Results of the edgeR Analysis."""
- section.markdown(logcpm_caption)
- section.plot(logcpm_vs_logfc)
+
+ dge = pd.read_csv(dge_file, sep="\t")
+ dge['sig'] = None
+ dge.loc[(dge["logFC"] > 0) & (dge['PValue'] < pval_thresh), 'sig'] = 'up'
+ dge.loc[(dge["logFC"] <= 0) & (dge['PValue'] < pval_thresh), 'sig'] = 'down'
+ dge.loc[(dge["PValue"] >= pval_thresh), 'sig'] = 'not_sig'
+
+ plot = scatterplot(
+ data=dge, x='logCPM', y='logFC', hue='sig',
+ palette=['#E32636', '#7E8896', '#0A22DE'],
+ hue_order=['up', 'not_sig', 'down'], marker='circle')
+ plot._fig.x_range.start = 10
+ plot._fig.xaxis.axis_label = "Average log CPM"
+ plot._fig.yaxis.axis_label = "Log-fold change"
+ plot.legend = dict(orient='horizontal', top=30)
+ # Should opacity of the symbols be lowered?
+ plot._fig.title = "Average copy per million (CPM) vs Log-fold change (LFC)"
+ EZChart(plot)
-def salmon_table(salmon_counts, section):
+def salmon_table(salmon_counts):
"""Create salmon counts summary table."""
salmon_counts = pd.read_csv(salmon_counts, sep='\t')
salmon_counts.set_index("Reference", drop=True, append=False, inplace=True)
salmon_size_top = salmon_counts.sum(axis=1).sort_values(ascending=False)
salmon_counts = salmon_counts.applymap(np.int64)
- salmon_count_caption = """
- Table showing the annotated Transcripts Per Million
+ h5("Transcripts Per Million")
+ p("""Table showing the annotated Transcripts Per Million
identified by Minimap2 mapping and Salmon transcript
- detection with the highest
- number of mapped reads"""
- section.markdown("### Transcripts Per Million ")
- section.markdown(salmon_count_caption, "salmon-head-caption")
- section.table(
- salmon_counts.loc[salmon_size_top.index].head(n=100), index=True)
+ detection. Displaying the top 100 transcripts with the highest
+ number of mapped reads""")
+
+ salmon_counts = salmon_counts[sorted(salmon_counts.columns)]
+ DataTable.from_pandas(
+ salmon_counts.loc[salmon_size_top.index].head(n=100), use_index=True)
def get_translations(gtf):
"""Create dict with gene_name and gene_references."""
- fn = open(gtf).readlines()
- gene_txid = {}
- gene_geid = {}
- geid_gname = {}
+ with open(gtf) as fh:
+ gene_txid = {}
+ gene_geid = {}
+ geid_gname = {}
- def get_feature(row, feature):
- return row.split(feature)[1].split(
- ";")[0].replace('=', '').replace("\"", "").strip()
+ def get_feature(row, feature):
+ return row.split(feature)[1].split(
+ ";")[0].replace('=', '').replace("\"", "").strip()
- for i in fn:
- if i.startswith("#"):
- continue
- # Different gtf/gff formats contain different attributes
- # and different formating (eg. gene_name="xyz" or gene_name "xyz")
- if 'gene_name' in i:
- gene_name = get_feature(i, "gene_name")
- elif 'gene_id' in i:
- gene_name = get_feature(i, 'gene_id')
- elif 'gene' in i:
- gene_name = get_feature(i, "gene")
- else:
- continue
+ for i in fh:
+ if i.startswith("#"):
+ continue
+ # Different gtf/gff formats contain different attributes
+ # and different formating (eg. gene_name="xyz" or gene_name "xyz")
+ gene_name = None
+ for var_name in ["gene_name", "gene_id", "gene"]:
+ if var_name in i:
+ gene_name = get_feature(i, var_name)
+ break
- if 'ref_gene_id' in i:
- gene_reference = get_feature(i, 'ref_gene_id')
- elif 'gene_id' in i:
- gene_reference = get_feature(i, 'gene_id')
- else:
- gene_reference = gene_name
- if 'transcript_id' in i:
- transcript_id = get_feature(i, 'transcript_id')
- else:
- transcript_id = "unknown"
- if 'gene_id' in i:
- gene_id = get_feature(i, 'gene_id')
- else:
- gene_id = gene_name
- gene_txid[transcript_id] = gene_name
- gene_geid[gene_id] = gene_reference
- geid_gname[gene_reference] = gene_name
+ if 'ref_gene_id' in i:
+ gene_reference = get_feature(i, 'ref_gene_id')
+ elif 'gene_id' in i:
+ gene_reference = get_feature(i, 'gene_id')
+ else:
+ gene_reference = gene_name
+ if 'transcript_id' in i:
+ transcript_id = get_feature(i, 'transcript_id')
+ else:
+ transcript_id = "unknown"
+ if 'gene_id' in i:
+ gene_id = get_feature(i, 'gene_id')
+ else:
+ gene_id = gene_name
+ gene_txid[transcript_id] = gene_name
+ gene_geid[gene_id] = gene_reference
+ geid_gname[gene_reference] = gene_name
return gene_txid, gene_geid, geid_gname
def de_section(
stringtie, dge, dexseq, dtu,
tpm, report, filtered, unfiltered,
- gene_counts):
+ gene_counts, aln_stats_dir, pval_threshold=0.01):
"""Differential expression sections."""
- section = report.add_section()
- section.markdown("# Differential expression.")
- section.markdown("""
-This section shows differential gene expression
-and differential isoform usage. Salmon was used to
-assign reads to individual annotated isoforms defined by
-the GTF-format annotation.
-These counts were used to perform a statistical analysis to identify
-the genes and isoforms that show differences in abundance between
-the experimental conditions.
-Any novel genes or transcripts that do not have relevant gene or transcript IDs
-are prefixed with MSTRG for use in differential expression analysis.
-Find the full sequences of any transcripts in the
-`final_non_redundant_transcriptome.fasta` file.
- """)
- section.markdown("### Alignment summary stats")
- alignment_stats = pool_csvs("seqkit")
- alignment_summary_df = create_summary_table(alignment_stats)
- alignment_summary_df = alignment_summary_df.fillna(0).applymap(np.int64)
- section.table(alignment_summary_df, key='alignment-stats', index=True)
- salmon_table(tpm, section)
- gene_txid, gene_name, geid_gname = get_translations(stringtie)
- # Use dictionaries to add gene names to the counts tsv files to help users
- copy_dge_add_gene_names(dge, geid_gname, "results_dge.tsv")
- copy_dge_add_gene_names(gene_counts, geid_gname, "all_gene_counts.tsv")
- copy_filtered_add_gene_names(
- filtered, geid_gname, "filtered_transcript_counts_with_genes.tsv")
- copy_filtered_add_gene_names(
- unfiltered, geid_gname, "unfiltered_transcript_counts_with_genes.tsv")
- copy_tpm_add_gene_names(tpm, gene_txid, "unfiltered_tpm_transcript_counts.tsv")
+ with report.add_section("Differential expression", "DE"):
- # Add tables to report
- dge_section(dge, section, gene_name)
- dexseq_section(dexseq, section, gene_name)
- dtu_section(dtu, section, gene_txid, gene_name)
- # missing dtu plots at the moment as too many
- section.markdown("""
-### View dtu_plots.pdf file to see plots of differential isoform usage
-""")
+ p("""This section shows differential gene expression
+ and differential isoform usage. Salmon was used to
+ assign reads to individual annotated isoforms defined by
+ the GTF-format annotation.
+ These counts were used to perform a statistical analysis to identify
+ the genes and isoforms that show differences in abundance between
+ the experimental conditions.
+ Any novel genes or transcripts that do not have relevant gene or transcript IDs
+ are prefixed with MSTRG for use in differential expression analysis.
+ 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)
+ h5("Alignment summary stats")
+ alignment_summary_df.index.name = "statistic"
+ DataTable.from_pandas(alignment_summary_df, use_index=True)
+
+ salmon_table(tpm)
+ gene_txid, gene_name, geid_gname = get_translations(stringtie)
+
+ # Add gene names columns to counts files and write out
+ # for publishing to user dir.
+ df_dge = pd.read_csv(dge, sep='\t')
+ df_dge.insert(0, 'gene_name', df_dge.index.map(lambda x: geid_gname.get(x)))
+ df_dge.to_csv('results_dge.tsv', index=True, index_label="gene_id", sep="\t")
+
+ # write_dge(gene_counts, geid_gname, "all_gene_counts.tsv")
+ df_gene_counts = pd.read_csv(gene_counts, sep='\t')
+ df_gene_counts.insert(
+ 0, 'gene_name', df_gene_counts.index.map(lambda x: geid_gname.get(x)))
+ df_gene_counts.to_csv(
+ 'results_dge.tsv', index=True, index_label="gene_id", sep="\t")
+
+ df_filtered = pd.read_csv(filtered, sep='\t')
+ df_filtered.insert(1, "gene_name", df_filtered.gene_id.map(
+ lambda x: geid_gname.get(x)))
+ df_filtered.to_csv(
+ 'filtered_transcript_counts_with_genes.tsv', index=False, sep='\t')
+
+ df_unfiltered = pd.read_csv(unfiltered, sep='\t')
+ df_unfiltered.insert(1, "gene_name", df_unfiltered.gene_id.map(
+ lambda x: geid_gname.get(x)))
+ df_unfiltered.to_csv(
+ 'unfiltered_transcript_counts_with_genes.tsv', index=False, sep='\t')
+
+ df_tpm = pd.read_csv(tpm, sep='\t')
+ df_tpm.insert(1, "gene_name", df_tpm.Reference.map(
+ lambda x: gene_txid.get(x)))
+ df_tpm.to_csv("unfiltered_tpm_transcript_counts.tsv", index=False, sep='\t')
+
+ # Add tables to report
+ dge_section(dge, gene_name, pval_threshold)
+ dexseq_section(dexseq, gene_name, pval_threshold)
+ dtu_section(dtu, gene_txid, gene_name)
diff --git a/bin/workflow_glue/generate_tracking_summary.py b/bin/workflow_glue/generate_tracking_summary.py
deleted file mode 100755
index 02a0eef..0000000
--- a/bin/workflow_glue/generate_tracking_summary.py
+++ /dev/null
@@ -1,66 +0,0 @@
-#!/usr/bin/env python
-"""Generate per-transcript class sumarrarry files from gffcompare."""
-
-import os
-import sys
-
-import pandas as pd
-
-from .util import wf_parser # noqa: ABS101
-
-
-def argparser():
- """Argument parser for entrypoint."""
- parser = wf_parser("generate_tracking_summary")
- parser.add_argument("--tracking", required=True, help="")
- parser.add_argument("--output_dir", required=True, help="")
- parser.add_argument("--annotation", required=False, default=None, help="")
-
- return parser
-
-
-def generate_tracking_summary(tracking_file, output_dir, annotations=None):
- """Write per transcript class gffcompare tracking files."""
- # results/gffcompare/gff_comparestringtie.tracking
- write_empty_tsvs = False
- tracking_headings = [
- "query_transfrag_id", "query_locus_id", "ref_gene_id",
- "class", "details"]
- nice_names = {
- '=': 'complete', 'c': 'contained', 'k': 'containment',
- 'm': 'retained', 'n': 'retained (partial)', 'j': 'multi',
- 'e': 'single', 'o': 'overlap', 's': 'opposite',
- 'x': 'exonic', 'i': 'intron', 'y': 'contains', 'p': 'runon',
- 'r': 'repeat', 'u': 'unknown'}
-
- if os.path.exists(annotations):
- tracking = pd.read_csv(
- tracking_file, sep="\t", names=tracking_headings[1:],
- index_col=0)
-
- d = pd.DataFrame(tracking['class'].value_counts()) \
- .reset_index().rename(columns={'index': 'class', 'class': 'count'})
- d['description'] = [nice_names[x] for x in d['class']]
-
- # write a separate table for each class
- for class_code, table in tracking.groupby('class'):
- if not write_empty_tsvs and table.empty:
- sys.stdout("Skipping: No transcripts found for: {}".format(
- class_code))
- continue
-
- path = tracking_file + ".{}.tsv".format(class_code)
- table.to_csv(path)
- else:
- sys.stdout(
- "Skipping classification summary as no annotation provided.")
-
-
-def main(args):
- """Run entry point."""
- assert os.path.isfile(args.tracking)
- assert os.path.isdir(args.output_dir)
- if args.annotation:
- os.path.isfile(args.annotation)
- generate_tracking_summary(
- args.tracking, output_dir=args.output_dir, annotations=args.annotation)
diff --git a/bin/workflow_glue/make_isoform_table.py b/bin/workflow_glue/make_isoform_table.py
deleted file mode 100755
index bb46f88..0000000
--- a/bin/workflow_glue/make_isoform_table.py
+++ /dev/null
@@ -1,68 +0,0 @@
-#!/usr/bin/env python
-"""Make report tables and data for plotting."""
-from pathlib import Path
-
-import numpy as np
-import pandas as pd
-from .util import wf_parser # noqa: ABS101
-
-
-def argparser():
- """Argument parser for entrypoint."""
- parser = wf_parser("Prepare report data")
- parser.add_argument(
- '--sample_id', help="Sample ID", required=True)
- parser.add_argument(
- '--gffcompare_dir',
- help="The gffcompare output directory",
- required=False,
- type=Path)
-
- return parser
-
-
-def make_isoform_table(gffcompare_dir, sample_id):
- """Make an isoform summary table."""
- try:
- tmap_file = next(gffcompare_dir.glob('*.tmap'))
- except StopIteration:
- raise ValueError("Cannot find .tmap file in {}".format(gffcompare_dir))
- dtypes = {
- 'ref_gene_id': str,
- 'ref_id': str,
- 'class_code': str,
- 'qry_id': str,
- 'num_exons': np.uint16,
- 'cov': np.uint32,
- 'len': np.uint32
- }
- df = pd.read_csv(
- tmap_file, sep='\t+',
- index_col=None,
- usecols=list(dtypes.keys()),
- dtype=dtypes)
-
- if len(df) == 0: # No transcripts. Write a header only result file
- df = pd.DataFrame(
- columns=list(dtypes.keys()) + ['sample_id', 'parent gene iso num'])
- df.to_csv(f'{sample_id}_transcripts_table.tsv', sep='\t', index=False)
- else:
- df = df.assign(sample_id=sample_id)
-
- # Make a column of number of isoforms in parent gene
- gb = df.groupby(['ref_gene_id']).count()
- gb.rename(columns={'ref_id': 'num_isoforms'}, inplace=True)
-
- df['parent gene iso num'] = df.apply(
- lambda x: gb.loc[(x.ref_gene_id), 'num_isoforms'], axis=1)
-
- # Unclassified transcripts should not be lumped together
- df.loc[df.class_code == 'u', 'parent gene iso num'] = None
-
- df.to_csv(f'{sample_id}_transcripts_table.tsv', sep='\t', index=False)
-
-
-def main(args):
- """Entry point."""
- if args.gffcompare_dir:
- make_isoform_table(args.gffcompare_dir, args.sample_id)
diff --git a/bin/workflow_glue/parse_gffcompare.py b/bin/workflow_glue/parse_gffcompare.py
new file mode 100755
index 0000000..20ed5b2
--- /dev/null
+++ b/bin/workflow_glue/parse_gffcompare.py
@@ -0,0 +1,262 @@
+#!/usr/bin/env python
+"""Make report tables and data for plotting."""
+import os
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+from .util import get_named_logger, wf_parser # noqa: ABS101
+
+
+def argparser():
+ """Argument parser for entrypoint."""
+ parser = wf_parser("Parse gffcompare")
+ parser.add_argument(
+ '--sample_id', help="Sample ID", required=True)
+ parser.add_argument(
+ '--gffcompare_dir',
+ help="The gffcompare output directory",
+ required=False,
+ type=Path)
+ parser.add_argument(
+ '--isoform_table_out',
+ help="Output path for per-isoform table",
+ type=Path)
+ parser.add_argument(
+ '--tracking',
+ help="gffcompare tracking file",
+ type=Path)
+ parser.add_argument(
+ "--annotation",
+ required=False,
+ default=None, help="Reference annotation GFF file")
+
+ return parser
+
+
+def _parse_stat_line(sl):
+ """Parse a stats line."""
+ res = {}
+ tmp = sl.split(':')[1].split('|')
+ res['sensitivity'] = float(tmp[0].strip())
+ res['precision'] = float(tmp[1].strip())
+ return res
+
+
+def _parse_matching_line(line):
+ """Parse a matching line."""
+ tmp = line.split(':')[1].strip()
+ return int(tmp)
+
+
+def _parse_mn_line(line):
+ """Parse a miss or novel line."""
+ res = {}
+ tmp = line.split(':')[1].strip()
+ tmp = tmp.split('/')
+ res['value'] = int(tmp[0])
+ tmp = tmp[1].split('(')
+ res['value_total'] = int(tmp[0].strip())
+ res['percent'] = float(tmp[1].split('%)')[0])
+ return res
+
+
+def _parse_total_line(line):
+ """Parse a total line."""
+ res = {}
+ tmp = line.split(':')[1].strip()
+ tmp = tmp.split('in')
+ res['transcripts'] = int(tmp[0].strip())
+ tmp = tmp[1].split('loci')
+ res['loci'] = int(tmp[0].strip())
+ tmp = int(tmp[1].split('(')[1].split(' ')[0])
+ res['me_transcripts'] = tmp
+ return res
+
+
+def parse_gffcmp_stats(gffcompare_stats, sample_id, outpath):
+ """Parse a gffcompare stats file.
+
+ Gffcompare stats file
+ :param gffcompare_stats: Path to the gffcompare stats file.
+ :returns: Return as tuple of dataframes containing:
+ perfromance statistics, match statistics, miss statistics,
+ novel statistics, total statistics.
+ :rtype: tuple
+ """
+ performance = []
+ missed = []
+ novel = []
+ total = []
+
+ with open(gffcompare_stats, 'r') as fh:
+ for line in fh:
+ line = line.strip()
+ if len(line) == 0:
+ continue
+
+ # Parse totals:
+ if line.startswith('# Query mRNAs'):
+ r = _parse_total_line(line)
+ total.append([r['loci'], 'loci', 'query'])
+ total.append([r['transcripts'], 'transcripts', 'query'])
+ total.append([r['me_transcripts'], 'multexonic', 'query'])
+ if line.startswith('# Reference mRNAs '):
+ r = _parse_total_line(line)
+ total.append([r['loci'], 'loci', 'reference'])
+ total.append([r['transcripts'], 'transcripts', 'reference'])
+ total.append([r['me_transcripts'], 'multexonic', 'reference'])
+
+ # Parse basic statistics:
+ if line.startswith('Base level'):
+ st = _parse_stat_line(line)
+ performance.append((st['sensitivity'], 'Sensitivity', 'Base'))
+ performance.append((st['precision'], 'Precision', 'Base'))
+ if line.startswith('Exon level'):
+ st = _parse_stat_line(line)
+ performance.append((st['sensitivity'], 'Sensitivity', 'Exon'))
+ performance.append((st['precision'], 'Precision', 'Exon'))
+ if line.startswith('Intron level'):
+ st = _parse_stat_line(line)
+ performance.append((st['sensitivity'], 'Sensitivity', 'Intron'))
+ performance.append((st['precision'], 'Precision', 'Intron'))
+ if line.startswith('Intron chain level'):
+ st = _parse_stat_line(line)
+ performance.append((st['sensitivity'], 'Sensitivity', 'Intron_chain'))
+ performance.append((st['precision'], 'Precision', 'Intron_chain'))
+ if line.startswith('Transcript level'):
+ st = _parse_stat_line(line)
+ performance.append((st['sensitivity'], 'Sensitivity', 'Transcript'))
+ performance.append((st['precision'], 'Precision', 'Transcript'))
+ if line.startswith('Locus level'):
+ st = _parse_stat_line(line)
+ performance.append((st['sensitivity'], 'Sensitivity', 'Locus'))
+ performance.append((st['precision'], 'Precision', 'Locus'))
+
+ # Parse missing statistics:
+ if line.startswith('Missed exons'):
+ r = _parse_mn_line(line)
+ missed.append((r['value'], 'Missed', 'Exons'))
+ missed.append((r['value_total'], 'total', 'Exons'))
+ missed.append((r['percent'], 'Percent', 'Exons'))
+ if line.startswith('Missed introns'):
+ r = _parse_mn_line(line)
+ missed.append((r['value'], 'Missed', 'Introns'))
+ missed.append((r['value_total'], 'total', 'Introns'))
+ missed.append((r['percent'], 'Percent', 'Introns'))
+ if line.startswith('Missed loci'):
+ r = _parse_mn_line(line)
+ missed.append((r['value'], 'Missed', 'Loci'))
+ missed.append((r['value_total'], 'total', 'Loci'))
+ missed.append((r['percent'], 'Percent', 'Loci'))
+
+ # Parse novel statistics:
+ if line.startswith('Novel exons'):
+ r = _parse_mn_line(line)
+ novel.append((r['value'], 'Novel', 'Exons'))
+ novel.append((r['value_total'], 'Total', 'Exons'))
+ novel.append((r['percent'], 'Percent_novel', 'Exons'))
+ if line.startswith('Novel introns'):
+ r = _parse_mn_line(line)
+ novel.append((r['value'], 'Novel', 'Introns'))
+ novel.append((r['value_total'], 'Total', 'Introns'))
+ novel.append((r['percent'], 'Percent_novel', 'Introns'))
+ if line.startswith('Novel loci'):
+ r = _parse_mn_line(line)
+ novel.append((r['value'], 'Novel', 'Loci'))
+ novel.append((r['value_total'], 'Total', 'Loci'))
+ novel.append((r['percent'], 'Percent_novel', 'Loci'))
+
+ def write_records(records, fn):
+ pd.DataFrame.from_records(records, columns=['counts', 'type', 'source']) \
+ .to_csv(outpath / fn, sep='\t')
+
+ write_records(total, 'Totals.tsv')
+ write_records(missed, 'Missed.tsv')
+ write_records(performance, 'Performance.tsv')
+ write_records(novel, 'Novel.tsv')
+
+
+def tracking_summary(tracking_file, output_dir, annotations=None):
+ """Write per transcript class gffcompare tracking files."""
+ tracking_headings = [
+ "query_transfrag_id", "query_locus_id", "ref_gene_id",
+ "class", "details"]
+ nice_names = {
+ '=': 'complete', 'c': 'contained', 'k': 'containment',
+ 'm': 'retained', 'n': 'retained (partial)', 'j': 'multi',
+ 'e': 'single', 'o': 'overlap', 's': 'opposite',
+ 'x': 'exonic', 'i': 'intron', 'y': 'contains', 'p': 'runon',
+ 'r': 'repeat', 'u': 'unknown'}
+
+ if os.path.exists(annotations):
+ tracking = pd.read_csv(
+ tracking_file, sep="\t", names=tracking_headings[1:],
+ index_col=0)
+
+ df = (
+ pd.DataFrame(tracking['class'].value_counts())
+ .reset_index()
+ .rename(columns={'index': 'class', 'class': 'Count'})
+ )
+
+ df['Percent'] = round(df['Count'] * 100 / df['Count'].sum(), 2)
+ df['description'] = [nice_names[x] for x in df['class']]
+
+ df = df.sort_values('Count', ascending=True)
+ df.to_csv(output_dir / 'tracking_summary.tsv', sep='\t')
+
+ else:
+ logger = get_named_logger('trackingSum')
+ logger.info("Skipping classification summary as no annotation provided.")
+
+
+def make_isoform_table(gffcompare_dir, sample_id, outpath):
+ """Make an isoform summary table."""
+ try:
+ tmap_file = next(gffcompare_dir.glob('*.tmap'))
+ except StopIteration:
+ raise ValueError("Cannot find .tmap file in {}".format(gffcompare_dir))
+ dtypes = {
+ 'ref_gene_id': str,
+ 'ref_id': str,
+ 'class_code': str,
+ 'qry_id': str,
+ 'num_exons': np.uint16,
+ 'cov': np.uint32,
+ 'len': np.uint32
+ }
+ df = pd.read_csv(
+ tmap_file, sep='\t+',
+ index_col=None,
+ usecols=list(dtypes.keys()),
+ dtype=dtypes)
+
+ if df.empty: # No transcripts. Write a header only result file
+ df = pd.DataFrame(
+ columns=list(dtypes.keys()) + ['sample_id', 'parent gene iso num'])
+ df.to_csv(f'{sample_id}_transcripts_table.tsv', sep='\t', index=False)
+ else:
+ df = df.assign(sample_id=sample_id)
+
+ # Make a column of number of isoforms in parent gene
+ gb = df.groupby(['ref_gene_id']).count()
+ gb.rename(columns={'ref_id': 'num_isoforms'}, inplace=True)
+
+ df['parent gene iso num'] = df.apply(
+ lambda x: gb.loc[(x.ref_gene_id), 'num_isoforms'], axis=1)
+
+ # Unclassified transcripts should not be lumped together
+ df.loc[df.class_code == 'u', 'parent gene iso num'] = None
+
+ df.to_csv(outpath, sep='\t', index=False)
+
+
+def main(args):
+ """Entry point."""
+ if args.gffcompare_dir: # TODO: should this every be optional?
+ stats = args.gffcompare_dir / 'str_merged.stats'
+ parse_gffcmp_stats(stats, args.sample_id, args.gffcompare_dir)
+ make_isoform_table(args.gffcompare_dir, args.sample_id, args.isoform_table_out)
+ tracking_summary(
+ args.tracking, args.gffcompare_dir, args.annotation)
diff --git a/bin/workflow_glue/reheader_samstream.py b/bin/workflow_glue/reheader_samstream.py
new file mode 100755
index 0000000..1d6528a
--- /dev/null
+++ b/bin/workflow_glue/reheader_samstream.py
@@ -0,0 +1,296 @@
+"""Reheader a SAM in a stream.
+
+When using the bam2fq -> minimap2 pattern for (re)aligning BAM data, we
+lose any existing RG and PG headers. This is particularly egregious when
+handling basecalled data as lines related to dorado basecalling settings
+as well as dorado RG headers are lost; orphaning RG tags in the reads.
+This is problematic for downstream anaylses that would like to read the
+XAM header to intelligently determine how to handle the reads based on
+the basecaller model and basecaller configuration.
+
+This script handles:
+ - Inserting RG, PG and CO lines from an existing XAM header into the
+ header of the SAM emitted from minimap2's alignment stream
+ - Inserting a PG header to indicate that a call to bam2fq was made
+ - Updating the first streamed PG.PP parent tag with the last PG.ID
+ of the existing XAM header to maintain a chain of custody
+ - Updating any streamed PG.ID (and PG.PP) tags to avoid collisions
+ with inserted PG.ID
+
+Handling collisions may seem like overkill but it is anticipated that
+this script will be called immediately after minimap2, any previous
+attempt to use minimap2 will lead to ambiguity. This would be the
+expected case where users have used wf-basecalling or wf-alignment to
+align a set of reads, only to realign them to another reference (eg.
+via wf-human-variation). Arguably, we should remove older references to
+minimap2 as they will have been invalidated by the call to bam2fq but
+removing PG records and sticking the PG chain back together seems more
+fraught with annoying future bugs than simply resolving conflicts.
+
+This script will explode on a stream that contains:
+ - PG lines in the original header where the last PG in the chain is
+ ambiguous, or where the parent PP IDs are not injective
+ - PG lines in the stream that do not appear in the order of their
+ chain (that is if a PG.PP refers to a PG.ID that has not been
+ encountered yet)
+
+SQ lines are retained after an HD line. That is to say, the most recent
+set of SQ lines observed after an HD will appear in the final output.
+SQ, RG, PG and CO lines are emitted as a group together, with elements
+written out in the order observed.
+
+PG lines are naively appended to the last PG element in the chain. No
+attempt is made to keep multiple program chains intact as this can lead
+to bloated headers. Broken PG metadata is a known problem (see
+samtools/hts-specs#275) but one that is preferable to headers that
+become unwieldly large to process: there IS an upper limit to a SAM
+header's size after all.
+
+This script takes advantage of minimap2's SAM output to immediately
+reheader the stream before any downstream calls to other programs pollute
+the PG header. This script is a little overkill but attempts to be robust
+with handling PG collisions and more obviously encapsulates reheadering
+behaviour, and leaves some room to do more clever things as necessary.
+"""
+import sys
+
+from .util import wf_parser # noqa: ABS101
+
+
+class SamHeader:
+ """An overkill container to manage merging PG lines in SAM headers.
+
+ Collision handling is simple. If a PG.ID is duplicated by the stream
+ then we add a suffix to its name and keep an eye out for the
+ corresponding PG.PP later. We assume that headers emitted by the
+ stream are chronological because this script should not be called as
+ part of any complicated pipework other than immediately following
+ minimap2.
+ """
+
+ def __init__(self):
+ """Initialise a collision aware PG container."""
+ self.remapped_pgids = {}
+ self.collision_suffix = 0
+
+ # Default HD, in case the new stream does not provide one
+ self.hd = "@HD\tVN:1.6\tSO:unknown"
+
+ # We'll merge RG, CO and PG
+ self.rg_records = []
+ self.co_records = []
+ self.pg_records = []
+
+ # We keep the most recently observed block of SQ records by
+ # resetting SQ on the first SQ seen after non-SQ. We cannot
+ # rely on HD being emitted (as minimap2 does not do this!)
+ self.sq_records = []
+ self.reset_sq = False
+
+ self.observed_rgids = set()
+ self.observed_pgids = set()
+ self.last_pgid = None
+
+ @staticmethod
+ def str_to_record(line):
+ """Return an appropriate struct for a given string record."""
+ try:
+ record_type, record_data = line.strip().split('\t', 1)
+ except ValueError:
+ raise Exception(f"Record type could not be determined: {line}")
+
+ if len(record_type) > 3:
+ raise Exception(f"Record type malformed: {record_type}")
+
+ record = {}
+ if record_type in ["@HD", "@CO", "@SQ"]:
+ return record_type, record_data
+ elif record_type in ["@RG", "@PG"]:
+ allowed_keys = {
+ "@RG": ["ID", "BC", "CN", "DS", "DT", "FO", "KS", "LB", "PG", "PI", "PL", "PM", "PU", "SM"], # noqa:E501
+ "@PG": ["ID", "PN", "CL", "PP", "DS", "VN"]
+ }
+ for field in record_data.strip().split('\t'):
+ k, v = field.split(':', 1)
+ if k not in allowed_keys[record_type]:
+ raise Exception(f"{record_type} with bad key '{k}': {record_data}")
+ record[k] = v
+ if "ID" not in record:
+ raise Exception(f"{record_type} with no ID: {record_data}")
+ return record_type, record
+ else:
+ raise Exception(f"Unknown record type: {line}")
+
+ @staticmethod
+ def record_to_str(record_type, record_data):
+ """Form a string from a header record."""
+ if record_type in ["@PG", "@RG"]:
+ tags = [f"{k}:{v}" for k, v in record_data.items()]
+ return f"{record_type}\t" + '\t'.join(tags)
+ elif record_type in ["@SQ", "@CO"]:
+ return f"{record_type}\t{record_data}"
+
+ @staticmethod
+ def resolve_pg_chain(pg_dicts):
+ """Check links between PG.ID and PP.ID, exploding if inconsistent."""
+ links = {}
+ # Document links between all ID and their PP parent
+ pgids_without_ppid = 0
+ for pgd in pg_dicts:
+ pgid = pgd["ID"]
+ pgpp = pgd.get("PP")
+ links[pgid] = pgpp
+ if pgpp is None:
+ pgids_without_ppid += 1
+ if len(links) > 0:
+ # If there are links, exactly one should have a None parent
+ # to indicate the first PG in the chain. Explode if we see
+ # no head or multiple heads.
+ if pgids_without_ppid == 0:
+ raise Exception("PG chain does not have a head.")
+ elif pgids_without_ppid > 1:
+ raise Exception("PG chain has multiple heads.")
+ for source in links:
+ head = source
+ path = [head]
+ while True:
+ head = links[head]
+ if head is None:
+ break
+ if head in path:
+ path.append(head)
+ raise Exception(f"PG chain appears to contain cycle: {path}")
+ path.append(head)
+ # This function is only really called to catch any explosions
+ # but we'll return the links here as it is useful for testing
+ return links
+
+ def _bump_pg_collider(self):
+ """Alter the collision suffix after determining a collision."""
+ self.collision_suffix += 1
+
+ def _uncollide_pgid(self, pgid):
+ """Return an uncollided string for a given PG ID."""
+ new_pgid = f"{pgid}-{self.collision_suffix}"
+ self.remapped_pgids[pgid] = new_pgid
+ self._bump_pg_collider()
+ return new_pgid
+
+ def add_line(self, line):
+ """Add a header line to the header."""
+ record_type, record = self.str_to_record(line)
+
+ if record_type == "@HD":
+ self.hd = f"@HD\t{record}"
+ elif record_type == "@CO":
+ self.co_records.append(record)
+ elif record_type == "@SQ":
+ if self.reset_sq:
+ self.sq_records = []
+ self.reset_sq = False
+ self.sq_records.append(record)
+ elif record_type == "@RG":
+ rgid = record["ID"]
+ if rgid not in self.observed_rgids:
+ self.observed_rgids.add(rgid)
+ self.rg_records.append(record)
+ elif record not in self.rg_records:
+ # if rgid has been seen before, abort if this record is different
+ raise Exception(
+ f"Duplicate RG with ID '{rgid}' conflicts with previously seen RG with same ID." # noqa:E501
+ )
+ elif record_type == "@PG":
+ pgid = record["ID"]
+ if pgid in self.observed_pgids:
+ # collision, rewrite the pgid
+ pgid = self._uncollide_pgid(pgid)
+ record["ID"] = pgid
+ else:
+ self.observed_pgids.add(pgid)
+
+ # maintain chain
+ ppid = record.get("PP")
+ if not ppid:
+ # record has no parent, this is either
+ # - the first record (last_pgid is None) so is the tail
+ # - an inserted record that needs its parent to be the current tail
+ if not self.last_pgid:
+ self.last_pgid = pgid
+ else:
+ record["PP"] = self.last_pgid
+ self.last_pgid = pgid
+ else:
+ if ppid not in self.observed_pgids:
+ raise Exception(
+ f"Encountered PG.PP '{ppid}' before observing corresponding PG.ID" # noqa:E501
+ )
+ # remap parent id (if needed)
+ record["PP"] = self.remapped_pgids.get(ppid, ppid)
+ # set tail to this record
+ self.last_pgid = pgid
+
+ self.pg_records.append(record)
+
+ if len(self.sq_records) > 0 and record_type != '@SQ':
+ self.reset_sq = True
+
+ return record
+
+ def write_header(self, fh):
+ """Write this header to a file handle."""
+ self.resolve_pg_chain(self.pg_records) # check PG header
+ fh.write(f"{self.hd}\n")
+ for sq in self.sq_records:
+ fh.write(self.record_to_str("@SQ", sq) + '\n')
+ for rg in self.rg_records:
+ fh.write(self.record_to_str("@RG", rg) + '\n')
+ for pg in self.pg_records:
+ fh.write(self.record_to_str("@PG", pg) + '\n')
+ for co in self.co_records:
+ fh.write(self.record_to_str("@CO", co) + '\n')
+
+
+def reheader_samstream(header_in, stream_in, stream_out, args):
+ """Run reheader_samstream."""
+ # read original header into container
+ sh = SamHeader()
+ for line in header_in:
+ sh.add_line(line)
+
+ # append user provided lines to container
+ for line in args.insert:
+ sh.add_line(line)
+
+ # read the header portion of the minimap2 stream
+ wrote_header = False
+ for line in stream_in:
+ if line[0] != '@':
+ # write out header on first alignment
+ sh.write_header(stream_out)
+ wrote_header = True
+ # and actually write the first alignment
+ stream_out.write(line)
+ break
+ sh.add_line(line)
+
+ # Pass through the rest of the alignments
+ for line in stream_in:
+ stream_out.write(line)
+ # If there were no alignments, we won't have hit the != @ case in the first stdin,
+ # and we won't have written the header out. Write a header if we haven't already.
+ if not wrote_header:
+ sh.write_header(stream_out)
+
+
+def argparser():
+ """Argument parser for entrypoint."""
+ parser = wf_parser("reheader_samstream")
+ parser.add_argument("header_in")
+ parser.add_argument("--insert", action="append", default=[])
+ return parser
+
+
+def main(args):
+ """reheader_samstream default entry point."""
+ with open(args.header_in) as header_in:
+ reheader_samstream(header_in, sys.stdin, sys.stdout, args)
diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py
index 0d55e6e..a09286b 100755
--- a/bin/workflow_glue/report.py
+++ b/bin/workflow_glue/report.py
@@ -1,28 +1,23 @@
#!/usr/bin/env python
"""Create workflow report."""
-from collections import Counter, defaultdict, OrderedDict
-import math
-import os
+import json
from pathlib import Path
-import sys
+import pickle
-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
-from aplanat.util import Colors
-from bokeh.layouts import gridplot
-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
-from bokeh.transform import dodge
-import gffutils
-import numpy as np
+from bokeh.models import HoverTool, Range1d
+from bokeh.models.tickers import AdaptiveTicker
+from dominate.tags import li, ul
+from dominate.util import raw
+import ezcharts as ezc
+from ezcharts.components import fastcat
+from ezcharts.components.ezchart import EZChart
+from ezcharts.components.reports.labs import LabsReport
+from ezcharts.components.theme import LAB_head_resources
+from ezcharts.layout.snippets import DataTable, Grid, Tabs
+from ezcharts.plots.categorical import barplot
+from ezcharts.util import get_named_logger
import pandas as pd
-import sigfig
-
from . import de_plots # noqa: ABS101
from .util import wf_parser # noqa: ABS101
@@ -31,7 +26,12 @@ def argparser():
"""Argument parser for entrypoint."""
parser = wf_parser("report")
parser.add_argument("--report", help="Report output file")
- parser.add_argument("--stats", help="Read stats files.", nargs='+')
+ parser.add_argument(
+ "--metadata", default='metadata.json', required=True,
+ help="sample metadata")
+ parser.add_argument(
+ "--stats", nargs='+',
+ help="Fastcat per-read stats, ordered as per entries in --metadata.")
parser.add_argument(
"--versions", required=True,
help="directory containing CSVs containing name,version.")
@@ -39,22 +39,19 @@ def argparser():
"--params", default=None, required=True,
help="A JSON file containing the workflow parameter key/values")
parser.add_argument(
- "--revision", default='unknown',
- help="git branch/tag of the executed workflow")
+ "--wf_version", default='unknown',
+ help="version of the executed workflow")
parser.add_argument(
- "--commit", default='unknown',
- help="git commit of the executed workflow")
- parser.add_argument(
- "--alignment_stats", required=False, default=None, nargs='*',
+ "--alignment_stats", required=False, default=None, type=Path,
help="TSV summary file of alignment statistics")
parser.add_argument(
- "--gff_annotation", required=False, nargs='+',
+ "--gff_annotation", required=False, type=Path,
help="transcriptome annotation gff file")
parser.add_argument(
- "--gffcompare_dir", required=False, default=None,
+ "--gffcompare_dir", required=False, default=None, type=Path,
help="gffcompare outout dir")
parser.add_argument(
- "--pychop_report", required=False, default=None,
+ "--pychop_report", required=False, default=None, type=Path,
help="TSV summary file of pychopper statistics")
parser.add_argument(
"--isoform_table", required=False, type=Path,
@@ -63,531 +60,149 @@ def argparser():
"--isoform_table_nrows", required=False, type=int, default=5000,
help="Maximum rows to display in isoforms table")
parser.add_argument(
- "--de_report", required=False, type=str, default=None,
+ "--transcriptome_summary", required=False, type=Path,
+ help="Path to dir containing transcriptome summary results ")
+ parser.add_argument(
+ "--de_report", required=False, type=Path, default=None,
help="Differential expression report optional")
parser.add_argument(
- "--de_stats", required=False, type=str, default=None, nargs='*',
+ "--de_stats", required=False, type=Path, default=None,
help="Differential expression report optional")
+ parser.add_argument(
+ "--pval_threshold", required=False, type=float, default=0.01,
+ help=(
+ "pvalue theshold for inclusion of differentially expressed genes"
+ " transcripts in plots "))
return parser
-def _parse_stat_line(sl):
- """Parse a stats line."""
- res = {}
- tmp = sl.split(':')[1]
- tmp = tmp.split('|')
- res['sensitivity'] = float(tmp[0].strip())
- res['precision'] = float(tmp[1].strip())
- return res
-
-
-def _parse_matching_line(line):
- """Parse a metching line."""
- tmp = line.split(':')[1].strip()
- return int(tmp)
-
-
-def _parse_mn_line(line):
- """Parse a miss or novel line."""
- res = {}
- tmp = line.split(':')[1].strip()
- tmp = tmp.split('/')
- res['value'] = int(tmp[0])
- tmp = tmp[1].split('(')
- res['value_total'] = int(tmp[0].strip())
- res['percent'] = float(tmp[1].split('%)')[0])
- return res
-
-
-def _parse_total_line(line):
- """Parse a total line."""
- res = {}
- tmp = line.split(':')[1].strip()
- tmp = tmp.split('in')
- res['transcripts'] = int(tmp[0].strip())
- tmp = tmp[1].split('loci')
- res['loci'] = int(tmp[0].strip())
- tmp = int(tmp[1].split('(')[1].split(' ')[0])
- res['me_transcripts'] = tmp
- return res
-
-
-def parse_gffcmp_stats(txt):
- """Parse a gffcompare stats file.
-
- :param txt: Path to the gffcompare stats file.
- :returns: Return as tuple of dataframes containing:
- perfromance statistics, match statistics, miss statistics,
- novel statistics, total statistics.
- :rtype: tuple
- """
- sensitivity = []
- precision = []
- level = []
-
- matching = OrderedDict()
-
- missed_level = []
- missed = []
- missed_total = []
- missed_percent = []
-
- novel_level = []
- novel = []
- novel_total = []
- novel_percent = []
-
- total_target = []
- total_loci = []
- total_transcripts = []
- total_multiexonic = []
-
- fh = open(txt, 'r')
- for line in fh:
- line = line.strip()
- if len(line) == 0:
- continue
- # Parse totals:
- if line.startswith('# Query mRNAs'):
- total_target.append('Query')
- r = _parse_total_line(line)
- total_loci.append(r['loci'])
- total_transcripts.append(r['transcripts'])
- total_multiexonic.append(r['me_transcripts'])
-
- if line.startswith('# Reference mRNAs '):
- total_target.append('Reference')
- r = _parse_total_line(line)
- total_loci.append(r['loci'])
- total_transcripts.append(r['transcripts'])
- total_multiexonic.append(r['me_transcripts'])
-
- # Parse basic statistics:
- if line.startswith('Base level'):
- st = _parse_stat_line(line)
- level.append('Base')
- sensitivity.append(st['sensitivity'])
- precision.append(st['precision'])
- if line.startswith('Exon level'):
- st = _parse_stat_line(line)
- level.append('Exon')
- sensitivity.append(st['sensitivity'])
- precision.append(st['precision'])
- if line.startswith('Intron level'):
- st = _parse_stat_line(line)
- level.append('Intron')
- sensitivity.append(st['sensitivity'])
- precision.append(st['precision'])
- if line.startswith('Intron chain level'):
- st = _parse_stat_line(line)
- level.append('Intron chain')
- sensitivity.append(st['sensitivity'])
- precision.append(st['precision'])
- if line.startswith('Transcript level'):
- st = _parse_stat_line(line)
- level.append('Transcript')
- sensitivity.append(st['sensitivity'])
- precision.append(st['precision'])
- if line.startswith('Locus level'):
- st = _parse_stat_line(line)
- level.append('Locus')
- sensitivity.append(st['sensitivity'])
- precision.append(st['precision'])
-
- # Parse match statistics:
- if line.startswith('Matching intron chains'):
- m = _parse_matching_line(line)
- matching['Intron chains'] = [m]
- if line.startswith('Matching transcripts'):
- m = _parse_matching_line(line)
- matching['Transcripts'] = [m]
- if line.startswith('Matching loci'):
- m = _parse_matching_line(line)
- matching['Loci'] = [m]
-
- # Parse missing statistics:
- if line.startswith('Missed exons'):
- missed_level.append('Exons')
- r = _parse_mn_line(line)
- missed.append(r['value'])
- missed_total.append(r['value_total'])
- missed_percent.append(r['percent'])
- if line.startswith('Missed introns'):
- missed_level.append('Introns')
- r = _parse_mn_line(line)
- missed.append(r['value'])
- missed_total.append(r['value_total'])
- missed_percent.append(r['percent'])
- if line.startswith('Missed loci'):
- missed_level.append('Loci')
- r = _parse_mn_line(line)
- missed.append(r['value'])
- missed_total.append(r['value_total'])
- missed_percent.append(r['percent'])
-
- # Parse novel statistics:
- if line.startswith('Novel exons'):
- novel_level.append('Exons')
- r = _parse_mn_line(line)
- novel.append(r['value'])
- novel_total.append(r['value_total'])
- novel_percent.append(r['percent'])
- if line.startswith('Novel introns'):
- novel_level.append('Introns')
- r = _parse_mn_line(line)
- novel.append(r['value'])
- novel_total.append(r['value_total'])
- novel_percent.append(r['percent'])
- if line.startswith('Novel loci'):
- novel_level.append('Loci')
- r = _parse_mn_line(line)
- novel.append(r['value'])
- novel_total.append(r['value_total'])
- novel_percent.append(r['percent'])
-
- fh.close()
-
- df_stats = pd.DataFrame(OrderedDict(
- [('Sensitivity', sensitivity), ('Precision', precision)]), index=level)
- df_match = pd.DataFrame(matching, index=['Matching'])
-
- df_miss = pd.DataFrame(
- OrderedDict(
- [('Total', missed_total),
- ('Missed', missed),
- ('Percent missed', missed_percent)]), index=missed_level)
-
- df_novel = pd.DataFrame(
- OrderedDict(
- [('Total', novel_total),
- ('Novel', novel),
- ('Percent novel', novel_percent)]), index=novel_level)
-
- df_total = pd.DataFrame(OrderedDict(
- [('Loci', total_loci), ('Transcripts', total_transcripts),
- ('Multiexonic', total_multiexonic)]), index=total_target)
-
- return df_stats, df_match, df_miss, df_novel, df_total
-
-
-def grouped_bar(df, title="", tilted_xlabs=False):
- """Create grouped bar plot from pandas dataframe.
-
- :param pandas.DataFrame
- Index:
- str: the x group labels - groups cluserted using these
- Columns:
- numeric: sub-groups of data - each sub group has same colour
-
- :returns bokaoh.plotting.figure instance
- """
- min_ = 0
- max_ = df.to_numpy().max()
- max_ = max_ + (max_ * 0.3) # Add some padding at top of plot for legends
- yrange = int(min_), int(max_)
-
- df['x_groups'] = df.index
- df = df.reset_index(drop=True)
- source = ColumnDataSource(data=df)
-
- p = figure(
- x_range=df['x_groups'], y_range=yrange, height=250, title=title,
- toolbar_location=None, tools="")
- i = 0
- # Use the dodge method to plot groups of bars
- # https://docs.bokeh.org/en/latest/docs/user_guide/categorical.html
- dodge_range = (-0.25, 0.25)
- current_dodge = dodge_range[0]
- 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
-
- for col in df.columns:
-
- num_colors = df.shape[1] - 1
- colors = list(zip(
- *[[Category10_10[x]] * (len(df.columns) - 1)
- for x in range(num_colors)]))
- colors = [item for sublist in colors for item in sublist]
-
- if col == 'x_groups':
- continue
- color = colors[i]
- i += 1
-
- width = df.size / 60
- v = p.vbar(
- x=dodge('x_groups', current_dodge, range=p.x_range), top=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 = "vertical"
- return p
-
-
def gff_compare_plots(report, gffcompare_outdirs):
"""Create various sections and plots in a WfReport.
- :param report: aplanat WFReport
:param gffcompare_outdirs: List of output directories from run_gffcompare
:return: None
- TODO: split this into separate functions
"""
- # If any of the gffcompare dirs are empty, skip this section
- if not all([any(Path(x).iterdir()) for x in gffcompare_outdirs]):
- return
-
# Plot overview panel:
- section = report.add_section()
- gffcompare_md = ('''
- ### Annotation summary
+ with report.add_section("Annotation summary", "Annotation"):
- The following plots summarize some of the output from
- [gffcompare](https://ccb.jhu.edu/software/stringtie/gffcompare.shtml)
+ raw("""The following plots summarize some of the output from
+ gffcompare
+ """)
- * **Totals**:
- Comparison of the number of stringtie-generated
- transcripts, multiexonic transcripts and loci between reference and query.
+ raw("""
+ - Totals: Comparison of the number of stringtie-generated
+ transcripts, multiexonic transcripts and loci between reference and
+ query
+ - Performance: How accurate are the query transcript annotations
+ with respect to the reference at various levels
+ - Missed: Features present in the reference, but absent in the
+ query
+ - Novel: Features present in the query transcripts,
+ but absent in the reference
""")
- * **Performance**:
- How accurate are the query transcript annotations with respect to the
- reference at various levels.
+ # Create the 4 bar plots of gffcompare summaries
+ gffcomp_fnames = ['Totals.tsv', 'Performance.tsv', 'Missed.tsv', 'Novel.tsv']
+ sample_tabs = Tabs()
- * **Missed**:
- Features present in the reference, but absent in the query
-
- * **Novel**:
- Features present in the query transcripts, but absent in the reference
- ''')
-
- tabs = []
- gff_fails = False
- sample_ids = []
- for dir_ in gffcompare_outdirs:
- sample_id = dir_.name
- sample_ids.append(sample_id) # Get sample ids fromt the folder name
- stats, _, miss, novel, total = \
- parse_gffcmp_stats(dir_ / 'str_merged.stats')
-
- if not any([x.empty for x in [stats, miss, novel, total]]):
- bar_totals = grouped_bar(total, title="Totals")
- bar_performance = grouped_bar(
- stats, title="Performance", tilted_xlabs=True)
- bar_missed = grouped_bar(miss, title="Missed")
- bar_novel = grouped_bar(novel, title="Novel")
- tabs.append(Panel(
- child=gridplot(
- [bar_totals, bar_performance, bar_missed, bar_novel],
- ncols=2, width=350, height=260), title=sample_id))
- else:
- 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:=',
- 'c': 'Contained:c',
- 'k': 'ReverseContained:k',
- 'm': 'RetainedIntron:m',
- 'n': 'PartRetainedIntron:n',
- 'j': 'PartialMatch:j',
- 'e': 'TransFragMatch:e',
- 's': 'OppositeMatch:s',
- 'o': 'OtherSameStrand:o',
- 'x': 'ExonicOpposite:o',
- 'y': 'RefInIntrons:y',
- 'p': 'PolymeraseRunon:p',
- 'r': 'Repeat:r',
- 'u': 'Intergenic:u',
- 'i': 'FullyIntronic:i',
- }
+ for dir_ in sorted(gffcompare_outdirs):
+ sample_id = dir_.name # Get sample ids from the folder name
+ with sample_tabs.add_tab(sample_id):
+ with Grid(columns=2):
+ for name in gffcomp_fnames:
+ df = pd.read_csv(dir_ / name, sep='\t')
+ if df.empty:
+ raw("""No gffcompare summary for this sample.
+ This could be due to incompatible reference fasta and gff
+ files""")
+ break
+ plot = barplot(data=df, x='source', y='counts', hue='type')
+ plot._fig.title = name.replace('.tsv', '')
+ EZChart(plot, height='250px')
# Plot overlaps panel:
- section = report.add_section()
- section.markdown('''
- ### Query transfrag classes
+ with report.add_section('Query transfrag classification', 'Classes'):
- The classes that are assigned by
- [gffcompare](https://ccb.jhu.edu/software/stringtie/gffcompare.shtml),
- which describe the relationship between query transfrag and the most
- similar reference transcript.
+ raw("""Summaries of the classes assigned by
+
+ gffcompare, which describe the relationship between query transfrag and
+ the most similar reference transcript.
- [This diagram](https://ccb.jhu.edu/software/stringtie/
- gffcompare_codes.png) illustrates the different classes.
- ''')
+
+ This diagram illustrates the different classes.
+ """)
+ tabs = Tabs()
+ for dir_ in sorted(gffcompare_outdirs):
+ sample_id = dir_.name # Get sample ids from the folder name
+ with tabs.add_tab(sample_id):
+ with Grid(columns=2):
+ track_file = dir_ / 'tracking_summary.tsv'
- tracking_dfs = []
+ df_tracking = (
+ pd.read_csv(track_file, sep="\t", index_col=0)
+ .sort_values('Count', ascending=True))
+ df_tracking.columns = [x.capitalize() for x in df_tracking.columns]
+ # Save class name for table
+ df_tracking.rename(columns={'Class': 'Class_tmp'}, inplace=True)
+ # Create new class name for barplot Description:Class
+ df_tracking['Class'] = \
+ df_tracking['Description'] + ":" + df_tracking['Class_tmp']
+ plot = barplot(
+ data=df_tracking, x='Count', y='Class', orient='h')
+ # revert to original class name
+ df_tracking = (
+ df_tracking
+ .drop(columns=['Class'])
+ .rename(columns={'Class_tmp': 'Class'}))
+ plot._fig.title.text = sample_id
+ plot._fig.yaxis.major_label_text_font_size = "10pt"
+ # Should the height of the plot be based on number of plot rows?
+ EZChart(plot, height='550px')
- track_files = [x / 'str_merged.tracking' for x in gffcompare_outdirs]
-
- df_tracking = load_data_add_sample_id(
- track_files, sample_ids,
- read_func=lambda x: pd.read_csv(
- x, sep="\t", header=None,
- usecols=[0, 3],
- names=['Count', 'Overlaps']))
-
- tabs = []
- for id_, df_track in df_tracking.groupby('sample_id'):
- tracking = df_track.groupby("Overlaps").count().reset_index()
- tracking.Overlaps = tracking.Overlaps.map(names)
- tracking['Percent'] = tracking.Count * 100 / tracking.Count.sum()
- tracking = tracking.sort_values("Count", ascending=False)
- track_bar = bars.simple_hbar(
- list(reversed(tracking['Overlaps'].values.tolist())),
- list(reversed(tracking['Percent'].values.tolist())),
- colors=Colors.cerulean, title=id_)
-
- tracking_dfs.append(tracking)
-
- tracking['Description'] = pd.Series(tracking.Overlaps.apply(
- lambda x: x.split(':')[0]))
-
- tracking['Code'] = pd.Series(tracking.Overlaps.apply(
- lambda x: x.split(':')[1]))
-
- tracking.drop(columns=['sample_id', 'Overlaps'], inplace=True)
- tracking = tracking[['Code', 'Description', 'Count', 'Percent']]
- tracking = tracking.round({
- 'Percent': 2
- })
-
- cols = [TableColumn(
- field=ci, title=ci, width=100) for ci in tracking.columns]
-
- track_table = DataTable(
- columns=cols, source=ColumnDataSource(tracking),
- index_position=None, width=500)
-
- tabs.append(Panel(
- child=gridplot([track_bar, track_table], ncols=2), title=id_)
- )
-
- cover_panel = Tabs(tabs=tabs)
- section.plot(cover_panel)
-
- def plot_isoforms_per_tpm_bin(
- df_code, class_code, sample_id, geomspace=False):
- """Make plots of number of isoforms per TPM coverage bin."""
- max_ = int(sigfig.round(df_code.TPM.max(), 2))
- if geomspace:
- bins = [math.ceil(x) for x in np.geomspace(10, max_, num=15)]
- else:
- bins = np.linspace(10, max_, 15)
- bins = np.unique(bins) # Low max_ can end up with duplicated bins
-
- groups = pd.cut(df_code.TPM, bins).value_counts()
- df_code.to_csv('dfcode.csv')
- df_temp = pd.DataFrame.from_dict(dict(
- x=[x.mid for x in groups.index], y=groups.values
- ))
- df_temp.sort_values(by='x', inplace=True)
-
- x = [str(math.ceil(x)) for x in df_temp.x]
- y = df_temp.y
-
- reads_per_iso_geom = \
- bars.simple_bar(x, y,
- title="{} - Num isoforms/TPM bin - "
- "gffcompare class code - '{}'".format(
- sample_id, class_code),
- colors=Colors.cerulean,
- x_axis_label='TPM',
- y_axis_label='Number of isoforms')
-
- reads_per_iso_geom.xaxis.major_label_orientation = math.pi / 2.8
- return reads_per_iso_geom
-
- log_plots = defaultdict(list)
-
- try:
- tmap_files = [next(x.glob('*.tmap')) for x in gffcompare_outdirs]
- except StopIteration:
- sys.stderr("Cannot find .tmap files in {}".format(gffcompare_outdirs))
- return
-
- df_tmap = load_data_add_sample_id(tmap_files, sample_ids)
-
- for id_, df in df_tmap.groupby('sample_id'):
-
- log_plots[id_].append(plot_isoforms_per_tpm_bin(
- df, 'all', id_, geomspace=True))
-
- for class_code, df_code in df.groupby('class_code'):
- log_plots[id_].append(plot_isoforms_per_tpm_bin(
- df_code, class_code, id_, geomspace=True))
-
- tabs = []
- for id_, sample_plots in log_plots.items():
- tabs.append(Panel(
- child=gridplot(sample_plots, ncols=2), title=id_))
-
- section.markdown('''
- ### Read coverage by gffcompare transfrag class''')
- cover_panel = Tabs(tabs=tabs)
- section.plot(cover_panel)
+ DataTable.from_pandas(
+ df_tracking.sort_values('Count', ascending=False),
+ use_index=False, paging=False, searchable=False)
def pychopper_plots(report, pychop_report):
"""Make plots from pychopper output.
- :param report: aplanat WFReport
+ :param report: ezcharts report
:param pychop_report: path to pychopper stats file
"""
- section = report.add_section()
- section.markdown('''
- ### Pychopper summary statisitcs
+ with report.add_section("Pychopper summary statisitcs", "Pychopper"):
+ raw("""The following plots summarize
+ pychopper output""")
- The following plots summarize [pychopper] output
- (https://github.com/epi2me-labs/pychopper)
+ ul(
+ li("""Full length: Reads with primers found in correct orientation at
+ both ends."""),
+ li("""Rescued: subreads extracted from chimeric reads"""),
+ li("""Unusable: Reads with missing or incorrect primer orientation"""),
+ li(""""+/-: Orientation of reads relative to the mRNA""")
+ )
- * **Pr.found**: Reads with primers found in correct orientation at
- both ends.
- * **Resc**: Reads 'rescued' from fused reads
- * **Unusable**: Read with missing or incorrect primer orientation
- * **+/-**: Orientation of reads relative to the mRNA
- ''')
+ df = pd.read_csv(pychop_report, sep='\t', index_col=0)
- plots = []
- df = pd.read_csv(pychop_report, sep='\t', index_col=0)
-
- for id_, df in df.groupby('sample_id'):
- df1 = df.set_index('Name', drop=True)
- df1 = df1.T[['Primers_found', 'Rescue', 'Unusable']]
- df1.rename(columns={
- 'Primers_found': 'Pr.found',
- 'Rescue': 'Resc',
- 'Unusable': 'Un'}, inplace=True)
- df2 = df[df.index == 'Strand']
- bar_chop = bars.simple_bar(
- df1.columns.values.tolist() + df2.Name.values.tolist(),
- df1.iloc[0].values.tolist() + df2.Value.values.tolist(),
- title='{} - Pychopper stats'.format(id_),
- colors=Colors.cerulean)
-
- plots.extend([bar_chop])
- grid = gridplot(
- plots, ncols=4, width=300, height=300)
- section.plot(grid)
+ with Grid(columns=2):
+ for sample, df in df.groupby('sample_id'):
+ df = df.set_index('Name', drop=True)
+ df = df.loc[['Primers_found', 'Rescue', 'Unusable', '+', '-']]
+ df = df.rename(
+ index={
+ 'Primers_found': 'full length',
+ 'Rescue': 'Rescued'
+ },
+ columns={
+ 'Value': 'n_reads'
+ }
+ ).reset_index(drop=False)
+ plot = barplot(
+ data=df, x='Name', y='n_reads')
+ plot._fig.title.text = f'{sample}'
+ plot._fig.xaxis.axis_label = ""
+ EZChart(plot, height="300px")
def transcript_table(report, isoform_table, max_rows):
@@ -595,182 +210,121 @@ def transcript_table(report, isoform_table, max_rows):
:param isoform_table: path to folder of isoform table files
"""
- section = report.add_section()
+ with report.add_section("Isoforms", "Isoforms"):
+ raw("""
+ This table details each isoforms identified per sample.
- # Should we 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
+ 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))
+ tabs = Tabs()
+ # drop some columns for the big table and do some filtering
- # drop some columns for the big table and do some filtering
- section.markdown('''
- ### Isoforms table
+ dfs = [pd.read_csv(x, sep='\t') for x in Path(isoform_table).iterdir()]
+ df = pd.concat(dfs)
+ for sample, df in df.groupby('sample_id'):
+ with tabs.add_tab(sample):
+ # 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, :]
- 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))
-
- dfs = []
- for file_ in Path(isoform_table).iterdir():
- d = pd.read_csv(file_, sep='\t')
- dfs.append(d)
- df = pd.concat(dfs)
-
- # 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)
+ # Sort by transcripts with the highest isoform diversity
+ df.sort_values('parent gene iso num', inplace=True, ascending=False)
+ DataTable.from_pandas(df, use_index=False)
-def transcriptome_summary(report, gffs):
+def transcriptome_summary(report, summaries_dir):
+ """Plot transcriptome summaries.
+
+ This section consists of four plots
+ 1: Isoforms per gene histogram
+ 2: Exons per transcript histogram
+ 3: Transcript lengths box plot
+ 4: Transcriptome summary table
"""
- Plot transcriptome summaries.
+ with report.add_section("Transcriptome summary", 'Summary'):
- Some of this data is available via gffcompare output, but the de novo
- pipeline skips that, so we do it all here.
+ plot_height = "300px"
+ tabs = Tabs()
+ data = {}
+ # Load all the dataframes upfront to get sample_id for sorting.
+ for summ_file in summaries_dir.glob('summary_*.tsv'):
+ with open(summ_file, 'rb') as fh:
+ summ = pickle.load(fh)
+ sample_id = summ['sample_id']
+ data[sample_id] = summ
- :param report: aplanat WFReport
- :param gffs: list of paths to gff transcriptome annotations
- """
- # test.db gets written to the git repo.
- section = report.add_section()
- section.markdown('''
- ### Transcriptome summary
- ''')
+ for sample_id in sorted(data):
+ summ = data[sample_id]
+ with tabs.add_tab(sample_id):
+ with Grid(columns=4):
- tabs = []
- for gff in gffs:
- sample_id = Path(gff).name
+ df_isoforms_per_gene = pd.DataFrame(
+ summ['isoforms_per_gene'].items(),
+ columns=['n_isoforms', 'count'])
+ iso_per_gene_plt = ezc.histplot(
+ data=df_isoforms_per_gene['n_isoforms'],
+ weights=df_isoforms_per_gene['count'],
+ discrete=df_isoforms_per_gene['n_isoforms'].max() < 20,
+ bins=20)
+ iso_per_gene_plt._fig.title = "Isoforms per gene"
+ iso_per_gene_plt._fig.xaxis.axis_label = 'Number of isoforms'
+ iso_per_gene_plt._fig.yaxis.axis_label = 'Number of genes'
+ iso_per_gene_plt._fig.xaxis.ticker = \
+ AdaptiveTicker(min_interval=1)
+ EZChart(iso_per_gene_plt, height=plot_height)
- plots = []
+ df_exons_per_transcript = pd.DataFrame(
+ summ['exons_per_transcript'].items(),
+ columns=['n_exons', 'count'])
+ exons_per_tr_plt = ezc.histplot(
+ data=df_exons_per_transcript['n_exons'],
+ weights=df_exons_per_transcript['count'],
+ discrete=df_exons_per_transcript['n_exons'].max() < 20,
+ bins=20)
+ exons_per_tr_plt._fig.title = "Exons per transcript"
+ exons_per_tr_plt._fig.xaxis.axis_label = 'Number of exons'
+ exons_per_tr_plt._fig.yaxis.axis_label = 'Number of transcripts'
+ exons_per_tr_plt._fig.xaxis.ticker = \
+ AdaptiveTicker(min_interval=1)
+ EZChart(exons_per_tr_plt, height=plot_height)
- db = gffutils.create_db(
- gff, dbfn=':memory:', force=True, keep_order=True,
- merge_strategy='merge', sort_attribute_values=True
- )
+ df_tr_len = pd.DataFrame.from_dict({
+ 'transcript_lengths': summ['transcript_lengths']})
+ df_tr_len['group'] = sample_id
+ tr_lens_plt = ezc.boxplot(
+ data=df_tr_len, x='group', y='transcript_lengths')
+ tr_lens_plt._fig.title = "Transcript lengths"
+ tr_lens_plt._fig.y_range = (
+ Range1d(0, df_tr_len['transcript_lengths'].max()))
- num_transcripts = db.count_features_of_type('transcript')
- num_genes = db.count_features_of_type('gene')
+ # Disable hover tool as ezcharts sometimes produces strange
+ # results for boxplots
+ hover_tools = tr_lens_plt._fig.select(dict(type=HoverTool))
+ for hover_tool in hover_tools:
+ tr_lens_plt._fig.tools.remove(hover_tool)
+ EZChart(tr_lens_plt, height=plot_height)
- transcript_lens = []
- exons_per_transcript = Counter()
- isoforms_per_gene = []
-
- for g in db.features_of_type('gene'):
-
- n_isos = len(list(db.children(g, featuretype='transcript')))
- isoforms_per_gene.append(n_isos)
-
- for t in db.children(
- g, featuretype='transcript', order_by='start'):
- tr_len = 0
- exons = list(enumerate(db.children(t, featuretype='exon')))
- if len(exons) == 0:
- continue
- for nx, ex in exons:
- tr_len += abs(ex.end - ex.start)
-
- exons_per_transcript[nx] += 1
-
- transcript_lens.append(tr_len)
-
- bar_isos = hist.histogram(
- [isoforms_per_gene], colors=[Colors.cerulean],
- title="isoforms per gene", binwidth=1)
- 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)
-
- box = bars.boxplot_series(
- [sample_id] * len(transcript_lens), transcript_lens,
- width=70, ylim=(min(transcript_lens), max(transcript_lens)),
- title='transcript lengths')
- plots.append(box)
-
- x, y = zip(*sorted(exons_per_transcript.items()))
-
- 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)
-
- df_sum = pd.DataFrame.from_dict(
- {'Total genes': [num_genes],
- 'Total transcripts': [num_transcripts],
- 'Max trans. len': max(transcript_lens),
- 'Min trans. len': min(transcript_lens)}).T
- df_sum.reset_index(drop=False, inplace=True)
-
- df_sum.columns = [' ', 'count']
- cols = [TableColumn(
- field=ci, title=ci, width=80) for ci in df_sum.columns]
- data_table = DataTable(
- columns=cols, source=ColumnDataSource(df_sum),
- index_position=None, width=180)
- plots.append(data_table)
-
- tabs.append(Panel(
- child=gridplot(plots, ncols=4,
- width=300, height=300), title=sample_id))
-
- cover_panel = Tabs(tabs=tabs)
- section.plot(cover_panel)
+ df_table = pd.DataFrame.from_dict(summ['summaries']).T
+ df_table = df_table.reset_index(drop=False)
+ df_table.columns = ['Transcriptome summary', '']
+ DataTable.from_pandas(
+ df_table, use_index=False, searchable=False, paging=False)
-def load_data_add_sample_id(files, sample_ids, read_func=None):
- """Load CSVs and concat into single into dataframe, and assign sample_id column."""
- df_ = pd.DataFrame()
- if not files:
- return None
- for id_, x in zip(sample_ids, files):
- if read_func:
- d = read_func(x)
- else:
- d = pd.read_csv(x, sep='\t+')
- d['sample_id'] = id_
- df_ = pd.concat([df_, d])
- return df_
-
-
-def seq_stats_tabs(report, stats):
- """Make tabs of sequence summaries by sample."""
- tabs = {}
- for summary_fn in stats:
- df_sample = pd.read_csv(summary_fn, sep="\t")
- sample_id = df_sample['sample_name'].iloc[0]
- rlp = read_length_plot(df_sample)
- rqp = read_quality_plot(df_sample)
- grid = gridplot(
- [rlp, rqp], ncols=2, sizing_mode="stretch_width")
- tabs[sample_id] = Panel(child=grid, title=sample_id)
- section = report.add_section()
- section.markdown("""
- ### Sequence summaries""")
- section.plot(Tabs(tabs=[tabs.get(x) for x in sorted(tabs)]))
-
-
-def de_section(report):
+def de_section(report, de_report_dir, de_aln_stats_dir, pval_threshold):
"""Make differential transcript expression section."""
- dexseq = os.path.join("de_report", "results_dexseq.tsv")
- dge = os.path.join("de_report", "results_dge.tsv")
- dtu = os.path.join("de_report", "results_dtu_stageR.tsv")
- stringtie = os.path.join("de_report", "stringtie_merged.gtf")
- tpm = os.path.join("de_report", "unfiltered_tpm_transcript_counts.tsv")
- filtered = os.path.join(
- "de_report", "filtered_transcript_counts_with_genes.tsv")
- unfiltered = os.path.join(
- "de_report", "unfiltered_transcript_counts_with_genes.tsv")
- gene_counts = os.path.join("de_report", "all_gene_counts.tsv")
+ dexseq = de_report_dir / "results_dexseq.tsv"
+ dge = de_report_dir / "results_dge.tsv"
+ dtu = de_report_dir / "results_dtu_stageR.tsv"
+ # GFF file can have gtf or gff extension
+ stringtie = next(de_report_dir.glob("*.g*f*"))
+ tpm = de_report_dir / "unfiltered_tpm_transcript_counts.tsv"
+ filtered = de_report_dir / "filtered_transcript_counts_with_genes.tsv"
+ unfiltered = de_report_dir / "unfiltered_transcript_counts_with_genes.tsv"
+ gene_counts = de_report_dir / "all_gene_counts.tsv"
# This will also add a gene name column to the above counts tsv files
de_plots.de_section(
stringtie=stringtie,
@@ -781,39 +335,57 @@ def de_section(report):
report=report,
filtered=filtered,
unfiltered=unfiltered,
- gene_counts=gene_counts)
+ gene_counts=gene_counts,
+ aln_stats_dir=de_aln_stats_dir,
+ pval_threshold=pval_threshold
+ )
def main(args):
"""Run the entry point."""
- report = WFReport(
- "Transcript isoform report", "wf-transcriptomes",
- revision=args.revision, commit=args.commit)
+ logger = get_named_logger("Report")
+ logger.info('Building report')
- seq_stats_tabs(report, args.stats)
+ report = LabsReport(
+ 'Workflow transcriptomes report', 'wf-transcriptomes',
+ args.params, args.versions, args.wf_version,
+ head_resources=[*LAB_head_resources])
+
+ with open(args.metadata) as metadata:
+ sample_details = [{
+ 'sample': d['alias'],
+ 'type': d['type'],
+ 'barcode': d['barcode']
+ } for d in json.load(metadata)]
+
+ with report.add_section('Read summary', 'Read summary'):
+ names = tuple(d['sample'] for d in sample_details)
+ stats = tuple(args.stats)
+ if len(stats) == 1:
+ stats = stats[0]
+ names = names[0]
+ fastcat.SeqSummary(stats, sample_names=names)
if args.alignment_stats is not None:
- stats_dfs = []
- for stats_file in args.alignment_stats:
- df = pd.read_csv(stats_file, sep='\t+')
- stats_dfs.append(df)
- aln_stats_df = pd.concat(stats_dfs)
- section = report.add_section()
- section.markdown('''
- ### Read mapping summary
+ with report.add_section('Alignment summary', 'Alignment'):
+ raw("""Reference genome alignment statistics. These were generated using
+ seqkit )
+ `seqkit bam -s` """)
- Summary of minimap2 mapping from
- [seqkit](https://bioinf.shenwei.me/seqkit/)
- `seqkit bam -s`''')
-
- section.table(aln_stats_df)
+ stats_dfs = []
+ for stats_file in args.alignment_stats.iterdir():
+ df = pd.read_csv(stats_file, sep='\\s+')
+ stats_dfs.append(df)
+ aln_stats_df = pd.concat(stats_dfs).sort_values(by='sample_id')
+ DataTable.from_pandas(aln_stats_df, use_index=False, export=True)
if args.pychop_report is not None:
- pychopper_plots(report, args.pychop_report)
+ # report is single file concatenated for all samples
+ pychopper_plots(report, next(args.pychop_report.iterdir()))
- # Results
+ # Do upstream
if args.gff_annotation is not None:
- transcriptome_summary(report, args.gff_annotation)
+ transcriptome_summary(report, args.transcriptome_summary)
if args.gffcompare_dir is not None:
gff_compare_plots(
@@ -824,12 +396,7 @@ def main(args):
transcript_table(report, args.isoform_table, args.isoform_table_nrows)
if args.de_report:
- de_section(report)
-
- # Arguments and software versions
- report.add_section(
- section=scomponents.version_table(args.versions))
- report.add_section(
- section=scomponents.params_table(args.params))
+ de_section(report, args.de_report, args.de_stats, pval_threshold=0.01)
report.write(args.report)
+ logger.info('Report writing finished')
diff --git a/bin/workflow_glue/summarise_gff.py b/bin/workflow_glue/summarise_gff.py
new file mode 100755
index 0000000..84483f9
--- /dev/null
+++ b/bin/workflow_glue/summarise_gff.py
@@ -0,0 +1,74 @@
+"""Get summary statistics from GFF file."""
+from collections import Counter
+from pathlib import Path
+import pickle
+
+import gffutils
+from .util import wf_parser # noqa: ABS101
+
+
+def argparser():
+ """Argument parser for entrypoint."""
+ parser = wf_parser("summ_gff")
+ parser.add_argument(
+ "gff",
+ help="Report output file",
+ type=Path)
+ parser.add_argument(
+ "sample_id",
+ help="Output TSV file path")
+ parser.add_argument(
+ "out",
+ default="gff_summary.tsv",
+ help="Output TSV file path",
+ type=Path)
+
+ return parser
+
+
+def main(args):
+ """Entry point."""
+ db = gffutils.create_db(
+ str(args.gff), dbfn=':memory:', force=True, keep_order=True,
+ merge_strategy='merge', sort_attribute_values=True
+ )
+
+ num_transcripts = db.count_features_of_type('transcript')
+ num_genes = db.count_features_of_type('gene')
+
+ transcript_lens = []
+ exons_per_transcript = Counter()
+ isoforms_per_gene = Counter()
+
+ for gene in db.features_of_type('gene'):
+
+ n_isos = len(list(db.children(gene, featuretype='transcript')))
+ isoforms_per_gene[n_isos] += 1
+
+ for transcript in db.children(
+ gene, featuretype='transcript', order_by='start'):
+ tr_len = 0
+ exons = list(db.children(transcript, featuretype='exon'))
+ if len(exons) == 0:
+ continue
+ exons_per_transcript[len(exons)] += 1
+ for ex in exons:
+ tr_len += abs(ex.end - ex.start)
+
+ transcript_lens.append(tr_len)
+
+ results = {
+ 'sample_id': args.sample_id,
+ 'summaries': {
+ 'Total genes': [num_genes],
+ 'Total transcripts': [num_transcripts],
+ 'Max trans. len': max(transcript_lens),
+ 'Min trans. len': min(transcript_lens)
+ },
+ 'transcript_lengths': transcript_lens,
+ 'exons_per_transcript': exons_per_transcript,
+ 'isoforms_per_gene': isoforms_per_gene
+ }
+
+ with open(args.out, 'wb') as fh:
+ pickle.dump(results, fh)
diff --git a/bin/workflow_glue/tests/test_test.py b/bin/workflow_glue/tests/test_test.py
deleted file mode 100755
index 496bc2a..0000000
--- a/bin/workflow_glue/tests/test_test.py
+++ /dev/null
@@ -1,10 +0,0 @@
-"""A dummy test."""
-
-import argparse
-
-from workflow_glue import report
-
-
-def test():
- """Just showing that we can import using the workflow-glue."""
- assert isinstance(report.argparser(), argparse.ArgumentParser)
diff --git a/docs/07_outputs.md b/docs/07_outputs.md
index 33ecfa0..9a46245 100644
--- a/docs/07_outputs.md
+++ b/docs/07_outputs.md
@@ -3,13 +3,13 @@ Output files may be aggregated including information for all samples or provided
| Title | File path | Description | Per sample or aggregated |
|-------|-----------|-------------|--------------------------|
| workflow report | wf-transcriptomes-report.html | a HTML report document detailing the primary findings of the workflow | aggregated |
-| Per file read stats | fastq_ingress_results/reads/fastcat_stats/per-file-stats.tsv | A TSV with per file read stats, including all samples. | aggregated |
-| Read stats | fastq_ingress_results/reads/fastcat_stats/per-read-stats.tsv | A TSV with per read stats, including all samples. | aggregated |
-| Run ID's | fastq_ingress_results/reads/fastcat_stats/run_ids | List of run IDs present in reads. | aggregated |
-| Meta map json | fastq_ingress_results/reads/metamap.json | Metadata used in workflow presented in a JSON. | aggregated |
-| Concatenated sequence data | fastq_ingress_results/reads/{{ alias }}.fastq.gz | Per sample reads concatenated in to one FASTQ file. | per-sample |
-| Assembled transcriptome | {{ alias }}_transcriptome.fas | Per sample assembled transcriptome. | per-sample |
-| Annotated assembled transcriptome | {{ alias }}_merged_transcriptome.fas | Per sample annotated assembled transcriptome. | per-sample |
+| Per file read stats | fastq_ingress_results/{{ alias }}//reads/fastcat_stats/per-file-stats.tsv | A TSV with per file read stats, including all samples. | aggregated |
+| Read stats | fastq_ingress_results/{{ alias }}//reads/fastcat_stats/per-read-stats.tsv | A TSV with per read stats, including all samples. | aggregated |
+| Run ID's | fastq_ingress_results/{{ alias }}//reads/fastcat_stats/run_ids | List of run IDs present in reads. | aggregated |
+| Meta map json | fastq_ingress_results/{{ alias }}//reads/metamap.json | Metadata used in workflow presented in a JSON. | aggregated |
+| Concatenated sequence data | fastq_ingress_results/{{ alias }}//reads/{{ alias }}.fastq.gz | Per sample reads concatenated in to one FASTQ file. | per-sample |
+| Assembled transcriptome | {{ alias }}_transcriptome.fas | Per sample assembled transcriptome. Not output if a reference annotation was supplied | per-sample |
+| Annotated assembled transcriptome | {{ alias }}_merged_transcriptome.fas | Per sample annotated assembled transcriptome. Only output if a reference annotation was supplied | per-sample |
| Alignment summary statistics | {{ alias }}_read_aln_stats.tsv | Per sample alignment summary statistics. | per-sample |
| GFF compare results. | {{ alias }}_gffcompare | All GFF compare output files. | per-sample |
| Differential gene expression results | de_analysis/results_dge.tsv | This is a gene-level result file that describes genes and their probability of showing differential expression between experimental conditions. | aggregated |
diff --git a/lib/ingress.nf b/lib/ingress.nf
index 2931357..72abfb5 100644
--- a/lib/ingress.nf
+++ b/lib/ingress.nf
@@ -329,15 +329,19 @@ def xam_ingress(Map arguments)
// - too many aligned files to safely and quickly merge (`samtools merge` opens
// all files at the same time and some machines might have low limits for
// open file descriptors)
- // * to_merge: flatMap > sort > group > merge
+ // * to_sortmerge: flatMap > sort > group > merge
+ // * to_merge: flatMap > group > merge
// - between 1 and `N_OPEN_FILES_LIMIT` aligned files
- no_files: n_files == 0
+ no_files: \
+ n_files == 0
indexed: \
n_files == 1 && (meta["is_unaligned"] || meta["is_sorted"]) && meta["src_xai"]
to_index: \
n_files == 1 && (meta["is_unaligned"] || meta["is_sorted"]) && !meta["src_xai"]
to_catsort: \
(n_files == 1) || (n_files > N_OPEN_FILES_LIMIT) || meta["is_unaligned"]
+ to_sortmerge: \
+ !meta["is_sorted"]
to_merge: true
}
@@ -345,6 +349,7 @@ def xam_ingress(Map arguments)
// only run samtools fastq on samples with at least one file
ch_to_fastq = ch_result.indexed.mix(
ch_result.to_index,
+ ch_result.to_sortmerge,
ch_result.to_merge,
ch_result.to_catsort
)
@@ -385,10 +390,13 @@ def xam_ingress(Map arguments)
}
// deal with samples with few-enough files for `samtools merge` first
- ch_merged = ch_result.to_merge
+ // we'll sort any unsorted files before merge
+ ch_merged = ch_result.to_sortmerge
| flatMap { meta, paths -> paths.collect { [meta, it] } }
| sortBam
+ | map { meta, bam, bai -> [meta, bam] } // drop index as merge does not need it
| groupTuple
+ | mix(ch_result.to_merge)
| mergeBams
| map{
meta, bam, bai ->
@@ -650,12 +658,13 @@ process validateIndex {
// Sort FOFN for samtools merge to ensure samtools sort breaks ties deterministically.
// Uses -c to ensure matching RG.IDs across multiple inputs are not unnecessarily modified to avoid collisions.
+// Note that samtools merge does not use the indexes so we do not provide them
process mergeBams {
label "ingress"
label "wf_common"
cpus 3
memory "4 GB"
- input: tuple val(meta), path("input_bams/reads*.bam"), path("input_bams/reads*.bam.bai")
+ input: tuple val(meta), path("input_bams/reads*.bam")
output: tuple val(meta), path("reads.bam"), path("reads.bam.bai")
script:
def merge_threads = Math.max(1, task.cpus - 1)
diff --git a/main.nf b/main.nf
index 8bd590d..ba2ac2c 100644
--- a/main.nf
+++ b/main.nf
@@ -264,7 +264,7 @@ process assemble_transcripts{
process merge_gff_bundles{
/*
- Merge gff bundles into a single gff file per sample.
+ Merge gff bundles into a single gff file per sample, and get summary statistics
*/
label 'isoforms'
cpus params.threads
@@ -274,6 +274,7 @@ process merge_gff_bundles{
tuple val(sample_id), path (gff_bundle)
output:
tuple val(sample_id), path("${sample_id}.gff"), emit: gff
+ tuple val(sample_id), path("transcriptome_summary.pickle"), emit: summary
script:
def merged_gff = "${sample_id}.gff"
"""
@@ -285,6 +286,11 @@ process merge_gff_bundles{
grep -v '#' \$fn >> $merged_gff
done
+
+ workflow-glue summarise_gff \
+ $merged_gff \
+ $sample_id \
+ transcriptome_summary.pickle
"""
}
@@ -315,17 +321,16 @@ process run_gffcompare{
gffcompare -o ${out_dir}/str_merged -r ${ref_annotation} \
${params.gffcompare_opts} ${query_annotation}
- workflow-glue generate_tracking_summary --tracking $out_dir/str_merged.tracking \
- --output_dir ${out_dir} --annotation ${ref_annotation}
-
mv *.tmap "${out_dir}"
mv *.refmap "${out_dir}"
cp "${out_dir}/str_merged.annotated.gtf" "${sample_id}_annotated.gtf"
- # Make an isoform table for report and user output.
- workflow-glue make_isoform_table \
+ workflow-glue parse_gffcompare \
--sample_id "${sample_id}" \
- --gffcompare_dir "${out_dir}"
+ --gffcompare_dir "${out_dir}" \
+ --isoform_table_out "${sample_id}_transcripts_table.tsv" \
+ --tracking $out_dir/str_merged.tracking \
+ --annotation ${ref_annotation}
"""
}
@@ -337,10 +342,9 @@ process get_transcriptome{
cpus 1
memory "2 GB"
input:
- tuple val(sample_id), path(transcripts_gff), path(gffcmp_dir), path(reference_seq)
+ tuple val(sample_id), path("transcripts.gff"), path(gffcompare_dir), path("reference.fa")
output:
- tuple val(sample_id), path("*.fas"), emit: transcriptome
-
+ tuple val(sample_id), path("*transcriptome.fas"), emit: transcriptome
script:
def transcriptome = "${sample_id}_transcriptome.fas"
def merged_transcriptome = "${sample_id}_merged_transcriptome.fas"
@@ -348,15 +352,11 @@ process get_transcriptome{
// so skip getting transcriptome FASTA from the annotated files.
if (params.ref_annotation){
"""
- gffread -g ${reference_seq} -w ${transcriptome} ${transcripts_gff}
- if [ "\$(ls -A $gffcmp_dir)" ];
- then
- gffread -F -g ${reference_seq} -w ${merged_transcriptome} $gffcmp_dir/str_merged.annotated.gtf
- fi
+ gffread -F -g reference.fa -w ${merged_transcriptome} $gffcompare_dir/str_merged.annotated.gtf
"""
} else {
"""
- gffread -g ${reference_seq} -w ${transcriptome} ${transcripts_gff}
+ gffread -g reference.fa -w ${transcriptome} "transcripts.gff"
"""
}
}
@@ -387,75 +387,59 @@ process merge_transcriptomes {
process makeReport {
- label "isoforms"
+ label "wf_common"
cpus 2
memory "4 GB"
-
+ publishDir "${params.out_dir}", mode: 'copy', pattern: "wf-transcriptomes-report.html"
input:
+ val metadata
+ path stats, stageAs: "stats_*"
path versions
+ val wf_version
path "params.json"
- path "pychopper_report/*"
- path "per_read_stats/?.gz"
- path "aln_stats/*"
- path "gffcmp_dir/*"
- path "gff_annotation/*"
- path "de_report/*"
- path "seqkit/*"
- path "isoforms_table/*"
+ path "transcriptome_aln_stats/*"
+ path pychopper, stageAs: "pychopper_report/*"
+ path aln_stats, stageAs: "aln_stats/*"
+ path gffcmp_dir, stageAs: "gffcmp_dir/*"
+ path gff_annotation, stageAs: "gff_annotation/*"
+ path de_report, stageAs: "de_report/*"
+ path isoforms_table, stageAs: "isoforms_table/*"
+ path "transcriptome_summary/summary_*.tsv"
+
output:
path ("wf-transcriptomes-*.html"), emit: report
- // If de analysis has been run output the counts files with gene name added.
path ("results_dge.tsv"), emit: results_dge, optional: true
path ("unfiltered_tpm_transcript_counts.tsv"), emit: tpm, optional: true
path ("unfiltered_transcript_counts_with_genes.tsv"), emit: unfiltered, optional: true
path ("filtered_transcript_counts_with_genes.tsv"), emit: filtered, optional: true
path ("all_gene_counts.tsv"), emit: gene_counts, optional: true
- shell:
- report_name = "wf-transcriptomes-report.html"
- '''
- if [ -f "de_report/OPTIONAL_FILE" ]; then
- dereport=""
- else
- dereport="--de_report true --de_stats "seqkit/*""
- mv de_report/*.g*f* de_report/stringtie_merged.gtf
- fi
- if [ -f "gff_annotation/OPTIONAL_FILE" ]; then
- OPT_GFF_ANNOTATION=""
- else
- OPT_GFF_ANNOTATION="--gff_annotation gff_annotation/*"
- fi
- if [ -f "gffcmp_dir/OPTIONAL_FILE" ]; then
- OPT_GFFCMP_DIR=""
- else
- OPT_GFFCMP_DIR="--gffcompare_dir gffcmp_dir/"
- fi
- if [ -f "aln_stats/OPTIONAL_FILE" ]; then
- OPT_ALN=""
- else
- OPT_ALN="--alignment_stats aln_stats/*"
- fi
- if [ -f "pychopper_report/OPTIONAL_FILE" ]; then
- OPT_PC_REPORT=""
- else
- OPT_PC_REPORT="--pychop_report pychopper_report/*"
- fi
- if [ -f "isoforms_table/OPTIONAL_FILE" ]; then
- OPT_ISO_TABLE=""
- else
- OPT_ISO_TABLE="--isoform_table isoforms_table"
- fi
- workflow-glue report --report !{report_name} \
- --versions !{versions} \
- --params params.json \
- ${OPT_ALN} \
- ${OPT_PC_REPORT} \
- --stats per_read_stats/* \
- ${OPT_GFF_ANNOTATION} \
- ${OPT_ISO_TABLE} \
- ${OPT_GFFCMP_DIR} \
- --isoform_table_nrows !{params.isoform_table_nrows} \
- ${dereport}
- '''
+ script:
+ String report_name = "wf-transcriptomes-report.html"
+ String metadata = new JsonBuilder(metadata).toPrettyString()
+ String gff_opts = gff_annotation.fileName.name == OPTIONAL_FILE.name ? "" : "--gff_annotation gff_annotation/"
+ String de_report_opts = de_report.fileName.name == OPTIONAL_FILE.name ? "" : "--de_report de_report/ --de_stats transcriptome_aln_stats/"
+ String gffcmp_opts = gffcmp_dir.fileName.name == OPTIONAL_FILE.name ? "" : "--gffcompare_dir gffcmp_dir/"
+ String aln_stats_opts = aln_stats.fileName.name == OPTIONAL_FILE.name ? "" : "--alignment_stats aln_stats/"
+ String pychop_opts = pychopper.fileName.name == OPTIONAL_FILE.name ? "" : "--pychop_report pychopper_report/"
+ String iso_table_opts = isoforms_table.fileName.name == OPTIONAL_FILE.name ? "" : "--isoform_table isoforms_table/"
+ """
+ echo '${metadata}' > metadata.json
+ workflow-glue report \
+ --report $report_name \
+ --versions $versions \
+ --wf_version $wf_version \
+ --params params.json \
+ $aln_stats_opts \
+ $pychop_opts \
+ --stats $stats \
+ --metadata metadata.json \
+ $gff_opts \
+ $iso_table_opts \
+ $gffcmp_opts \
+ --isoform_table_nrows ${params.isoform_table_nrows} \
+ $de_report_opts \
+ --transcriptome_summary transcriptome_summary/
+ """
}
@@ -598,9 +582,15 @@ workflow pipeline {
}
fastq_ingress_results = reads
- // replace `null` with path to optional file
- | map { [ it[0], it[1] ?: OPTIONAL_FILE, it[2] ?: OPTIONAL_FILE ] }
| collectFastqIngressResultsInDir
+ // fastq_ingress doesn't have the index; add one extra null for compatibility.
+ // We do not use variable name as assigning variable name with a tuple
+ // not matching (e.g. meta, bam, bai, stats <- [meta, bam, stats]) causes
+ // the workflow to crash.
+ reads = reads
+ .map{
+ it.size() == 4 ? it : [it[0], it[1], null, it[2]]
+ }
map_sample_ids_cls = {it ->
/* Harmonize tuples
output:
@@ -629,9 +619,8 @@ workflow pipeline {
String publish_bams = "BAMS"
software_versions = getVersions()
workflow_params = getParams()
- input_reads = reads.map{ meta, samples, stats -> [meta, samples]}
+ input_reads = reads.map{ meta, samples, index, stats -> [meta, samples]}
sample_ids = input_reads.flatMap({meta,samples -> meta.alias})
- per_read_stats = reads.map{ meta, samples, stats -> stats.resolve("per-read-stats.tsv.gz") }.toList()
if (!params.direct_rna){
preprocess_reads(input_reads)
@@ -682,7 +671,7 @@ workflow pipeline {
merge_gff = merge_gff_bundles.out.gff.map{ it -> it[1]}.collect()
// Output BAMS in a dedicated directory
- bam_results = assembly.bam.map{
+ bam_results = assembly.bam.map{
sample_id, bam, bai -> [bam, bai]}.flatten().map{ [it, publish_bams] }
}
else{
@@ -691,7 +680,6 @@ workflow pipeline {
merge_gff = OPTIONAL_FILE
assembly_stats = OPTIONAL_FILE
use_ref_ann = false
-
}
if (params.de_analysis){
sample_sheet = file(params.sample_sheet, type:"file")
@@ -716,24 +704,35 @@ workflow pipeline {
}
de = differential_expression(transcriptome, input_reads, sample_sheet, gtf)
de_report = de.all_de
- count_transcripts_file = de.count_transcripts
de_outputs = de.de_outputs
+ count_transcripts_file = de.count_transcripts
} else{
de_report = OPTIONAL_FILE
count_transcripts_file = OPTIONAL_FILE
}
-
+
+ // get metadata and stats files, keeping them ordered (could do with transpose I suppose)
+ reads.multiMap{ meta, path, index, stats ->
+ meta: meta
+ stats: stats
+ }.set { for_report }
+ metadata = for_report.meta.collect()
+ stats = for_report.stats.collect()
+
makeReport(
+ metadata,
+ stats,
software_versions,
+ workflow.manifest.version,
workflow_params,
+ count_transcripts_file,
pychopper_report,
- per_read_stats,
assembly_stats,
gff_compare,
merge_gff,
de_report,
- count_transcripts_file,
- isoforms_table)
+ isoforms_table,
+ merge_gff_bundles.out.summary.map {it[1]}.collect())
report = makeReport.out.report
@@ -769,6 +768,7 @@ workflow pipeline {
results = results.concat(de_results.map{ [it, "de_analysis"] })
}
+
results.concat(workflow_params.map{ [it, null]})
// IGV config
if (params.transcriptome_source == "precomputed" && params.igv){
@@ -822,7 +822,7 @@ workflow pipeline {
// get list of file names
// Absolute paths required for directories
igv_files = reads
- | map { meta, sample, stats -> meta.alias }
+ | map { meta, sample, index, stats -> meta.alias }
| toSortedList
| map { list -> list.collect{
[
@@ -830,9 +830,9 @@ workflow pipeline {
"$publish_bams/${it}_reads_aln_sorted.bam.bai"
]
} }
- | concat ( igv_index)
- | flatten
| concat (igv_ref)
+ | flatten
+ | concat ( igv_index)
| concat (gz_igv)
| flatten
| collectFile(name: "file-names.txt", newLine: true, sort: false)
@@ -844,7 +844,7 @@ workflow pipeline {
[displayMode: "SQUISHED", colorBy: "strand"], // bam extra opts
Channel.of(null), // vcf extra opts
)
-
+
results = results.concat(igv_conf.map{ [it, null]})
results = results.concat(bam_results)
}
diff --git a/output_definition.json b/output_definition.json
index c223e4e..1aedcb3 100644
--- a/output_definition.json
+++ b/output_definition.json
@@ -9,7 +9,7 @@
"type": "aggregated"
},
"read-stats-per-file": {
- "filepath": "fastq_ingress_results/reads/fastcat_stats/per-file-stats.tsv",
+ "filepath": "fastq_ingress_results/{{ alias }}//reads/fastcat_stats/per-file-stats.tsv",
"title": "Per file read stats",
"description": "A TSV with per file read stats, including all samples.",
"mime-type": "text/tab-separated-values",
@@ -17,7 +17,7 @@
"type": "aggregated"
},
"read-stats-per-read": {
- "filepath": "fastq_ingress_results/reads/fastcat_stats/per-read-stats.tsv",
+ "filepath": "fastq_ingress_results/{{ alias }}//reads/fastcat_stats/per-read-stats.tsv",
"title": "Read stats",
"description": "A TSV with per read stats, including all samples.",
"mime-type": "text/tab-separated-values",
@@ -25,7 +25,7 @@
"type": "aggregated"
},
"run-ids": {
- "filepath": "fastq_ingress_results/reads/fastcat_stats/run_ids",
+ "filepath": "fastq_ingress_results/{{ alias }}//reads/fastcat_stats/run_ids",
"title": "Run ID's",
"description": "List of run IDs present in reads.",
"mime-type": "text/txt",
@@ -33,7 +33,7 @@
"type": "aggregated"
},
"metamap": {
- "filepath": "fastq_ingress_results/reads/metamap.json",
+ "filepath": "fastq_ingress_results/{{ alias }}//reads/metamap.json",
"title": "Meta map json",
"description": "Metadata used in workflow presented in a JSON.",
"mime-type": "text/json",
@@ -41,7 +41,7 @@
"type": "aggregated"
},
"sample-data": {
- "filepath": "fastq_ingress_results/reads/{{ alias }}.fastq.gz",
+ "filepath": "fastq_ingress_results/{{ alias }}//reads/{{ alias }}.fastq.gz",
"title": "Concatenated sequence data",
"description": "Per sample reads concatenated in to one FASTQ file.",
"mime-type": "text/json",
@@ -51,7 +51,7 @@
"transcriptome": {
"filepath": "{{ alias }}_transcriptome.fas",
"title": "Assembled transcriptome",
- "description": "Per sample assembled transcriptome.",
+ "description": "Per sample assembled transcriptome. Not output if a reference annotation was supplied",
"mime-type": "text/x-fasta",
"optional": true,
"type": "per-sample"
@@ -59,7 +59,7 @@
"merged_transcriptome": {
"filepath": "{{ alias }}_merged_transcriptome.fas",
"title": "Annotated assembled transcriptome",
- "description": "Per sample annotated assembled transcriptome.",
+ "description": "Per sample annotated assembled transcriptome. Only output if a reference annotation was supplied",
"mime-type": "text/x-fasta",
"optional": true,
"type": "per-sample"
diff --git a/subworkflows/reference_assembly.nf b/subworkflows/reference_assembly.nf
index 412759b..734da48 100644
--- a/subworkflows/reference_assembly.nf
+++ b/subworkflows/reference_assembly.nf
@@ -28,11 +28,14 @@ process map_reads{
| samtools view -q ${params.minimum_mapping_quality} -F 2304 -Sb -\
| seqkit bam -j 1 -x -T '${ContextFilter}' -\
| samtools sort --write-index -@ 1 -o "${sample_id}_reads_aln_sorted.bam##idx##${sample_id}_reads_aln_sorted.bam.bai" - ;
- ((cat "${sample_id}_reads_aln_sorted.bam" | seqkit bam -s -j 1 - 2>&1) | tee ${sample_id}_read_aln_stats.tsv ) || true
+ ((cat "${sample_id}_reads_aln_sorted.bam" | seqkit bam -s -j 1 - 2>&1) | tee "${sample_id}_read_aln_stats.tsv" ) || true
- # Add sample id header and column
- sed "s/\$/${sample_id}/" "${sample_id}_read_aln_stats.tsv" \
- | sed "1 s/${sample_id}/sample_id/" > tmp
+ # Add sample id header and column; remove last column (File)
+ cat "${sample_id}_read_aln_stats.tsv" \
+ | sed "s/^/${sample_id} /" \
+ | sed "1 s/^${sample_id}/sample_id/" \
+ | awk 'NF{NF-=1};1' \
+ > tmp
mv tmp "${sample_id}_read_aln_stats.tsv"
if [[ -s "internal_priming_fail.tsv" ]];
