kfukuta/wf-transcriptomes-v202
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'bokehreport_CW-3077' into 'dev'

Browse Source 
Update report CW-3077

See merge request epi2melabs/workflows/wf-transcriptomes!168
This commit is contained in:
Neil Horner 2024-10-25 15:18:05 +00:00
commit 45a717aaf3
17 changed files with 1544 additions and 2083 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
CHANGELOG.md
View File

@ -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.

14
README.md
View File

@ -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 |

649
bin/workflow_glue/compute_cluster_quality.py
View File

@ -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()

454
bin/workflow_glue/configure_igv.py
View File

@ -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

513
bin/workflow_glue/de_plots.py
View File

@ -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)

66
bin/workflow_glue/generate_tracking_summary.py
View File

@ -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)

68
bin/workflow_glue/make_isoform_table.py
View File

@ -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)

262
bin/workflow_glue/parse_gffcompare.py Executable file
View File

@ -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)

296
bin/workflow_glue/reheader_samstream.py Executable file
View File

@ -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)

981
bin/workflow_glue/report.py
View File

File diff suppressed because it is too large Load Diff

74
bin/workflow_glue/summarise_gff.py Executable file
View File

@ -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)

10
bin/workflow_glue/tests/test_test.py
View File

@ -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)

14
docs/07_outputs.md
View File

@ -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 |

17
lib/ingress.nf
View File

@ -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)

170
main.nf
View File

@ -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)
@ -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)

14
output_definition.json
View File

@ -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"

11
subworkflows/reference_assembly.nf
View File

@ -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" ]];