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

Re-add ability to do direct RNA

Browse Source
This commit is contained in:
Neil Horner 2022-04-05 18:04:16 +00:00
parent 7dfea725a0
commit c2a0f8dc21
11 changed files with 202 additions and 204 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CHANGELOG.md
View File

@ -4,6 +4,10 @@ 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.0.0/),
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
# Unreleased
## Added
- --direct_rna option.
## [v0.1.1]
## Fixed
- Incorrect numbers and of transcripts caused by merging gff files with same gene and transcript ids

23
README.md
View File

@ -8,9 +8,9 @@ It has been adapted from two existing Snakemake pipelines:
* https://github.com/nanoporetech/pipeline-nanopore-denovo-isoforms
---
## Overview
* cDNA or direct RNA reads are initially and optionally preprocessed by [pychopper](https://github.com/nanoporetech/pychopper)
for identification of full-length reads, as well as trimming and orientation correction. <br>
* cDNA reads are initially preprocessed by [pychopper](https://github.com/nanoporetech/pychopper)
for identification of full-length reads, as well as trimming and orientation correction (This step is omited for
direct RNA reads)
Reference-based approach
* Full length reads are mapped to a supplied reference genome using [minimap2](https://github.com/lh3/minimap2) <br>
@ -22,11 +22,12 @@ using [gffcompare](http://ccb.jhu.edu/software/stringtie/gffcompare.shtml)
de novo-based approach (experimental!)
* Sequence clusters are generated using [isONclust2](https://github.com/nanoporetech/isONclust2)
* If a reference genome is supplied, cluster quality metrics are determined by comparing
with clusters generated from a minimap2 alignment
with clusters generated from a minimap2 alignment.
* A consensus sequence for each cluster is generated using [spoa](https://github.com/rvaser/spoa)
* Three rounds of polishing using racon and minimap2 to give a final polished CDS for each gene.
* Full-length reads are then mapped to these polished CDS
* Transcripts are assembled by stringtie as for the reference-based approach
* Full-length reads are then mapped to these polished CDS.
* Transcripts are assembled by stringtie as for the reference-based approach.
* Note: This is currently not supported with direct RNA reads.
For both approaches:
@ -67,7 +68,7 @@ to see the options for the workflow.
- Optional reference annotation in GFF2/3 format
**Example execution of a workflow for reference-based transcript assembly**
This uses a synthetic SIRV dataset so we need to tell minimap2 about the non-canonical spplice junctions with
This uses a synthetic SIRV dataset, so we need to tell minimap2 about the non-canonical splice junctions with
--minimap2_opts '-uf --splice-flank=no'
```
OUTPUT=~/output;
@ -94,6 +95,8 @@ or at the command line:<br>
- Threshold for including isoforms into interactive table `transcript_table_cov_thresh = 50`
- Run the denovo pipeline `denovo = true` (default false)
- To run the workflow with direct RNA reads `--direct_rna` (skips the pychopper step).
Pychopper and minimap2 can take options via `minimap2_opts` and `pychopper_opts`
<br>
@ -107,7 +110,7 @@ For example:
- pychopper can use one of two available backends for identifying primers in the raw reads
- nhmmscan `pychopper opts = '-m phmm'`
- edlib `pychopper opts = '-m edlib'`
<br>Note: edlib is used in the configs as it's quite a lot faster. However it may be less sensitive than nhmmscan.
<br>Note: edlib is set by default in the config as it's quite a lot faster. However it may be less sensitive than nhmmscan.
---
## Workflow outputs
@ -129,9 +132,7 @@ Each sample will also have it's own directory containing the following (dependen
* {sample_id}_transcriptome.fas
- A transcriptome derived from the query reads
* {sample_id}_merged_transcriptome.fas
- A transcriptome derived from the query reads + reference annotation
* merged_transcriptome.fas
-A transcriptome made from the combined query reads and reference annotation
- A transcriptome derived from the query reads + reference annotation
## Useful links

3
bin/check_sample_sheet.py
View File

@ -21,8 +21,7 @@ def main():
"Warning: sample sheet contains both 'alias' and "
'sample_id, using the former.')
samples['sample_id'] = samples['alias']
if 'barcode' not in samples.columns \
or 'sample_id' not in samples.columns:
if not set(['sample_id', 'barcode']).intersection(samples.columns):
raise IOError()
except Exception:
raise IOError(

59
bin/compute_cluster_quality.py
View File

@ -15,8 +15,9 @@ 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 sklearn.metrics.cluster import (
adjusted_rand_score, completeness_score,
homogeneity_score, v_measure_score)
matplotlib.use('Agg')
@ -368,14 +369,26 @@ def get_cluster_information(clusters, classes):
print("MIXED:", "Tot classes containing both:", len(
set(clustered_classes.keys()) & set(not_clustered_classes.keys())))
print("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)
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):
@ -400,14 +413,16 @@ def main(args):
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)
(
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")
@ -519,10 +534,10 @@ def main(args):
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')
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')

4
bin/generate_tracking_summary.py
View File

@ -60,8 +60,8 @@ def main(args):
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)
generate_tracking_summary(
args.tracking, output_dir=args.output_dir, annotations=args.annotation)
if __name__ == '__main__':

200
bin/report.py
View File

@ -23,52 +23,6 @@ import pandas as pd
import sigfig
def _vbar(x, top, title, **kwargs):
"""Vertical bar chart."""
fig = figure(title=title)
fig.vbar(x, top=top, **kwargs)
return fig
def _hbar(y, right, title='', fig_height=300, fig_width=300,
bar_height=0.1, **kwargs):
"""Horizontal bar chart."""
fig = figure(title=title, height=fig_height, width=fig_width)
yn = list(range(len(y)))
fig.hbar(yn,
right=right,
height=bar_height,
**kwargs)
# Overide the numerical labels with cate
fig.yaxis.ticker = yn
mapper = {k: v for (k, v) in zip(yn, y)}
fig.yaxis.major_label_overrides = mapper
return fig
def simple_hbar(df, y, right, title="", color=Colors.cerulean,
fig_kwargs={}, plot_kwargs={}):
"""Create a simple barplot.
:param groups: the grouping variable (the x-axis values).
:param values: the data for bars are drawn (the y-axis values).
:param kwargs: kwargs for bokeh figure.
"""
defaults = {
'output_backend': 'webgl',
'height': 300, 'width': 600}
defaults.update(fig_kwargs)
p = figure(y_range=df[y], height=250, title=title,
toolbar_location=None, tools="")
plot_kwargs.update({'height': 0.2})
p.hbar(y=df[y], right=df[right], **plot_kwargs)
return p
def _parse_stat_line(sl):
"""Parse a stats line."""
res = {}
@ -288,8 +242,9 @@ def grouped_bar(df, title="", tilted_xlabs=False):
df = df.reset_index(drop=True)
source = ColumnDataSource(data=df)
p = figure(x_range=df['x_groups'], y_range=yrange, height=250, title=title,
toolbar_location=None, tools="")
p = figure(
x_range=df['x_groups'], y_range=yrange, height=250, title=title,
toolbar_location=None, tools="")
i = 0
# Use the dodge method to plot groups of bars
# https://docs.bokeh.org/en/latest/docs/user_guide/categorical.html
@ -304,8 +259,9 @@ def grouped_bar(df, title="", tilted_xlabs=False):
for col in df.columns:
num_colors = df.shape[1] - 1
colors = list(zip(*[[Category10_10[x]] * (len(df.columns) - 1)
for x in range(num_colors)]))
colors = list(zip(
*[[Category10_10[x]] * (len(df.columns) - 1)
for x in range(num_colors)]))
colors = [item for sublist in colors for item in sublist]
if col == 'x_groups':
@ -314,9 +270,9 @@ def grouped_bar(df, title="", tilted_xlabs=False):
i += 1
width = df.size / 60
p.vbar(x=dodge('x_groups', current_dodge, range=p.x_range), top=col,
width=width, source=source, color=color,
legend_label=col)
p.vbar(
x=dodge('x_groups', current_dodge, range=p.x_range), top=col,
width=width, source=source, color=color, legend_label=col)
current_dodge += dodge_increment
p.x_range.range_padding = 0.1
@ -341,8 +297,7 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
# Plot overview panel:
section = report.add_section()
# TODO: update this based on current version
section.markdown('''
gffcompare_md = ('''
### Annotation summary
The following plots summarize some of the output from
@ -370,19 +325,26 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
stats, _, miss, novel, total = \
parse_gffcmp_stats(dir_ / 'str_merged.stats')
bar_totals = grouped_bar(total, title="Totals")
bar_performance = grouped_bar(stats, title="Performance",
tilted_xlabs=True)
bar_missed = grouped_bar(miss, title="Missed")
bar_novel = grouped_bar(novel, title="Novel")
if not any([x.empty for x in [stats, miss, novel, total]]):
bar_totals = grouped_bar(total, title="Totals")
bar_performance = grouped_bar(
stats, title="Performance", tilted_xlabs=True)
bar_missed = grouped_bar(miss, title="Missed")
bar_novel = grouped_bar(novel, title="Novel")
tabs.append(Panel(
child=gridplot(
[bar_totals, bar_performance, bar_missed, bar_novel],
ncols=2, width=350, height=260), title=id_))
tabs.append(Panel(
child=gridplot(
[bar_totals, bar_performance, bar_missed, bar_novel],
ncols=2, width=350, height=260), title=id_))
cover_panel = Tabs(tabs=tabs)
section.plot(cover_panel)
cover_panel = Tabs(tabs=tabs)
section.markdown(gffcompare_md)
section.plot(cover_panel)
else:
gffcompare_md = ('''
__Warning__: Gffcompare summary cannot be shown.
This could be due to incompatible reference fasta and gff
files.''')
section.markdown(gffcompare_md)
names = {
'=': 'ExactMatch:=',
@ -414,18 +376,19 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
[This diagram](https://ccb.jhu.edu/software/stringtie/
gffcompare_codes.png) illustrates the different classes.
''')
''')
tracking_dfs = []
print(gffcompare_outdirs)
track_files = [x / 'str_merged.tracking' for x in gffcompare_outdirs]
df_tracking = load_sample_data(track_files, sample_ids,
read_func=lambda x:
pd.read_csv(x, sep="\t", header=None,
usecols=[0, 3],
names=['Count', 'Overlaps']))
df_tracking = load_sample_data(
track_files, sample_ids,
read_func=lambda x: pd.read_csv(
x, sep="\t", header=None,
usecols=[0, 3],
names=['Count', 'Overlaps']))
tabs = []
for id_, df_track in df_tracking.groupby('sample_id'):
@ -433,9 +396,10 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
tracking.Overlaps = tracking.Overlaps.map(names)
tracking['Percent'] = tracking.Count * 100 / tracking.Count.sum()
tracking = tracking.sort_values("Overlaps")
track_bar = _hbar(
track_bar = bars.simple_hbar(
tracking['Overlaps'].values.tolist(),
tracking['Percent'].values.tolist(), title="{}".format(id_))
tracking['Percent'].values.tolist(),
colors=Colors.cerulean, title=id_)
tracking_dfs.append(tracking)
@ -448,13 +412,12 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
tracking.drop(columns=['sample_id', 'Overlaps'], inplace=True)
tracking = tracking[['Code', 'Description', 'Count', 'Percent']]
cols = [TableColumn(field=Ci, title=Ci, width=100)
for Ci in tracking.columns]
cols = [TableColumn(
field=Ci, title=Ci, width=100) for Ci in tracking.columns]
track_table = DataTable(columns=cols,
source=ColumnDataSource(tracking),
index_position=None,
width=500)
track_table = DataTable(
columns=cols, source=ColumnDataSource(tracking),
index_position=None, width=500)
tabs.append(Panel(
child=gridplot([track_bar, track_table], ncols=2), title=id_)
)
@ -462,13 +425,12 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
cover_panel = Tabs(tabs=tabs)
section.plot(cover_panel)
def plot_isoforms_per_tpm_bin(df_code, class_code, sample_id,
geomspace=False):
def plot_isoforms_per_tpm_bin(
df_code, class_code, sample_id, geomspace=False):
"""Make plots of number of isoforms per TPM coverage bin."""
max_ = int(sigfig.round(df_code.TPM.max(), 2))
if geomspace:
bins = [math.ceil(x) for x in
np.geomspace(10, max_, num=15)]
bins = [math.ceil(x) for x in np.geomspace(10, max_, num=15)]
else:
bins = np.linspace(10, max_, 15)
bins = np.unique(bins) # Low max_ can end up with duplicated bins
@ -527,11 +489,11 @@ def gff_compare_plots(report, gffcompare_outdirs: Path, sample_ids):
return df_tmap
def pychopper_plots(report, df):
def pychopper_plots(report, pychop_report):
"""Make plots from pychopper.cdna_classifier.py.
:param report: aplanat WFReport
:param df: DataFrame of Pychopper stats
:param pychop_report: path to pychopper stats file
"""
section = report.add_section()
section.markdown('''
@ -548,12 +510,15 @@ def pychopper_plots(report, df):
''')
plots = []
df = pd.read_csv(pychop_report, sep='\t', index_col=0)
for id_, df in df.groupby('sample_id'):
df1 = df.set_index('Name', drop=True)
df1 = df1.T[['Primers_found', 'Rescue', 'Unusable']]
df1.rename(columns={'Primers_found': 'Pr.found',
'Rescue': 'Resc',
'Unsable': 'Un'}, inplace=True)
df1.rename(columns={
'Primers_found': 'Pr.found',
'Rescue': 'Resc',
'Unusable': 'Un'}, inplace=True)
df2 = df[df.index == 'Strand']
bar_chop = bars.simple_bar(
df1.columns.values.tolist() + df2.Name.values.tolist(),
@ -562,8 +527,8 @@ def pychopper_plots(report, df):
colors=Colors.cerulean)
plots.extend([bar_chop])
grid = gridplot(plots, ncols=4,
width=300, height=300)
grid = gridplot(
plots, ncols=4, width=300, height=300)
section.plot(grid)
@ -602,8 +567,7 @@ def cluster_quality(cluster_qc_dir, report, sample_ids):
tabs = []
for id_, cluster_dir in zip(sample_ids, cluster_qc_dir):
plots = []
for fn in [
'v_ari_com_hom.csv', 'sing_nonsing.csv']:
for fn in ['v_ari_com_hom.csv', 'sing_nonsing.csv']:
# Skip the next two plots for now
# 'class_sizes1.csv', 'class_sizes2.csv']:
df = pd.read_csv(Path(cluster_dir) / fn)
@ -629,8 +593,9 @@ def transcript_table(report, df_tmaps, covr_threshold):
# all in single table and sample_id column? Currently it's the latter
# drop some columns for the big table and do some filtering
df = df_tmaps.drop(columns=['FPKM', 'qry_gene_id', 'major_iso_id',
'ref_match_len', 'TPM'])
df = df_tmaps.drop(
columns=[
'FPKM', 'qry_gene_id', 'major_iso_id', 'ref_match_len', 'TPM'])
df.sort_values('cov', ascending=True, inplace=True)
counts = list(range(len(df)))
@ -672,7 +637,7 @@ def transcript_table(report, df_tmaps, covr_threshold):
section.table(df, index=False)
def tanscriptome_summary(report, gffs, sample_ids, denovo=False):
def transcriptome_summary(report, gffs, sample_ids, denovo=False):
"""
Plot transcriptome summaries.
@ -745,12 +710,12 @@ def tanscriptome_summary(report, gffs, sample_ids, denovo=False):
if not denovo:
x, y = zip(*sorted(exons_per_transcript.items()))
ept_bar = _vbar(x, list(y),
title="Exons per transcript",
color=Colors.cerulean)
fig = figure(title="Exons per transcript")
fig.vbar(
x, top=list(y), color=Colors.cerulean)
ept_bar.xaxis.major_label_orientation = math.pi / 2.8
plots.append(ept_bar)
fig.xaxis.major_label_orientation = math.pi / 2.8
plots.append(fig)
df_sum = pd.DataFrame.from_dict(
{'Total genes': [num_genes],
@ -760,12 +725,11 @@ def tanscriptome_summary(report, gffs, sample_ids, denovo=False):
df_sum.reset_index(drop=False, inplace=True)
df_sum.columns = [' ', 'count']
cols = [TableColumn(field=Ci, title=Ci, width=80)
for Ci in df_sum.columns]
data_table = DataTable(columns=cols,
source=ColumnDataSource(df_sum),
index_position=None,
width=180)
cols = [TableColumn(
field=Ci, title=Ci, width=80) for Ci in df_sum.columns]
data_table = DataTable(
columns=cols, source=ColumnDataSource(df_sum),
index_position=None, width=180)
plots.append(data_table)
tabs.append(Panel(
@ -808,7 +772,6 @@ def main():
parser.add_argument(
"--commit", default='unknown',
help="git commit of the executed workflow")
parser.add_argument(
"--alignment_stats", required=False, default=None, nargs='*',
help="TSV summary file of alignment statistics")
@ -819,7 +782,7 @@ def main():
"--gffcompare_dir", required=False, default=None, nargs='*',
help="gffcompare outout dir")
parser.add_argument(
"--pychop_report", required=True, nargs='+',
"--pychop_report", required=False, default=None,
help="TSV summary file of pychopper statistics")
parser.add_argument(
"--sample_ids", required=True, nargs='+',
@ -860,8 +823,8 @@ def main():
section.table(df_aln_stats)
# workflow-specific plotting
tanscriptome_summary(report, args.gff_annotation, sample_ids,
denovo=args.denovo)
transcriptome_summary(
report, args.gff_annotation, sample_ids, denovo=args.denovo)
df_tmaps = gff_compare_plots(
report,
@ -870,17 +833,8 @@ def main():
report.write(args.report)
pc_df = pd.DataFrame()
for id_, pyc in zip(sample_ids, args.pychop_report):
try:
p = pd.read_csv(pyc, sep='\t', index_col=0)
except pd.errors.EmptyDataError:
continue
p['sample_id'] = id_
pc_df = pd.concat([pc_df, p])
if len(pc_df) > 0:
pychopper_plots(report, pc_df)
if args.pychop_report is not None:
pychopper_plots(report, args.pychop_report)
if df_tmaps is not None:
transcript_table(report, df_tmaps, args.transcript_table_cov_thresh)

33
bin/run_isonclust2.py
View File

@ -24,8 +24,8 @@ class Node:
def __repr__(self):
"""Get string repr of a node."""
return "Node:{} Level: {} File: {} Done: {} Left: {} Right: " \
"{} Parent: {}".format(
return "Node:{} Level: {} File: {} Done: " \
"{} Left: {} Right: {} Parent: {}".format(
self.Id, self.Level,
self.File, self.Done, self.Left.Id if
self.Left is not None else None,
@ -47,8 +47,10 @@ def build_job_tree():
"""Build a job tree of nodes."""
JOB_TREE = OrderedDict()
batches = glob("batches/isONbatch_*.cer")
batch_ids = [int(re.search('batches/isONbatch_(.*)\\.cer$', x).group(1))
for x in batches]
batch_ids = [
int(re.search(
'batches/isONbatch_(.*)\\.cer$', x).group(1))
for x in batches]
LEVELS = OrderedDict()
LEVELS[0] = []
for Id, bf in sorted(zip(batch_ids, batches), key=lambda x: x[0]):
@ -93,12 +95,14 @@ def main():
Path('clusters').mkdir(exist_ok=True)
job_tree, levels = build_job_tree()
init_template = 'isONclust2 cluster -x {} -v -Q -l batches/' \
'isONbatch_{}.cer -o clusters/isONcluster_{}.cer {}; ' \
'sync;\n'
template = 'isONclust2 cluster -x {} -v -Q -l clusters/isONcluster_{}' \
'.cer -r clusters/isONcluster_{}.cer -o clusters/isONcluster' \
'_{}.cer {}; sync\n'
init_template = (
'isONclust2 cluster -x {} -v -Q -l batches/isONbatch_{}.cer '
'-o clusters/isONcluster_{}.cer {}; '
'sync;\n')
template = (
'isONclust2 cluster -x {} -v -Q -l clusters/isONcluster_{}.cer '
'-r clusters/isONcluster_{}.cer -o clusters/isONcluster_{}.cer '
'{}; sync\n')
for nr, l in levels.items():
jobs_out = 'jobs_level_{}.sh'.format(nr)
@ -109,14 +113,15 @@ def main():
jr = init_template.format('sahlin', n.Id, n.Id, purge)
fh.write(jr)
else:
jr = template.format('sahlin', n.Left.Id,
n.Right.Id, n.Id, purge)
jr = template.format(
'sahlin', n.Left.Id, n.Right.Id, n.Id, purge)
fh.write(jr)
# Run a level in parallel
cmd = "parallel < {}".format(jobs_out)
sub.call(cmd, shell=True)
sub.call("ln -s `realpath clusters/isONcluster_{}"
".cer` isONcluster_ROOT.cer".format(n.Id), shell=True)
sub.call((
"ln -s `realpath clusters/isONcluster_{}.cer` "
"isONcluster_ROOT.cer".format(n.Id)), shell=True)
if __name__ == '__main__':

2
environment.yaml
View File

@ -6,7 +6,7 @@ channels:
- defaults
dependencies:
- python==3.8.*
- aplanat >=0.6.2
- aplanat >=0.6.4
- epi2melabs
- minimap2 ==2.24
- samtools ==1.14

70
main.nf
View File

@ -84,12 +84,16 @@ process preprocess_reads {
tuple val(sample_id), path(input_reads)
output:
tuple val(sample_id), path("${sample_id}_full_length_reads.fq"), emit: full_len_reads
tuple val(sample_id), path('*.tsv'), emit: report
path '*.tsv', emit: report
script:
"""
cdna_classifier.py -t ${params.threads} ${params.pychopper_opts} ${input_reads} ${sample_id}_full_length_reads.fq
mv cdna_classifier_report.tsv ${sample_id}_cdna_classifier_report.tsv
generate_pychopper_stats.py --data ${sample_id}_cdna_classifier_report.tsv --output .
# Add sample id column
sed "1s/\$/\tsample_id/; 1 ! s/\$/\t${sample_id}/" ${sample_id}_cdna_classifier_report.tsv > tmp
mv tmp ${sample_id}_cdna_classifier_report.tsv
"""
}
@ -266,12 +270,12 @@ process makeReport {
path versions
path "params.json"
val denovo
path pychopper_report
tuple val(sample_ids),
path(seq_summaries),
path(aln_stats),
path(gffcmp_dir),
path(cdna_class_report),
path(gff_annotation)
path(seq_summaries),
path(aln_stats),
path(gffcmp_dir),
path(gff_annotation)
output:
path("wf-isoforms-*.html"), emit: report
script:
@ -280,12 +284,13 @@ process makeReport {
def report_name = "wf-isoforms-report.html"
def OPT_ALN = denovo ? '' : "--alignment_stats ${aln_stats}"
def OPT_DENOVO = denovo ? "--denovo" : ''
def OPT_PC_REPORT = pychopper_report.name.startsWith('OPTIONAL_FILE') ? '' : "--pychop_report ${pychopper_report}"
"""
report.py --report $report_name \
--versions $versions \
--params params.json \
$OPT_ALN \
--pychop_report $cdna_class_report \
$OPT_PC_REPORT \
--sample_ids $sids \
--summaries $seq_summaries \
--gffcompare_dir $gffcmp_dir \
@ -345,16 +350,25 @@ workflow pipeline {
software_versions = getVersions()
workflow_params = getParams()
preprocess_reads(summariseConcatReads.out.input_reads)
if (!params.direct_rna){
preprocess_reads(summariseConcatReads.out.input_reads)
full_len_reads = preprocess_reads.out.full_len_reads
pychopper_report = preprocess_reads.out.report.collectFile(keepHeader: true)
}
else{
full_len_reads = summariseConcatReads.out.input_reads
pychopper_report = file("$projectDir/data/OPTIONAL_FILE")
}
if (params.denovo){
println("Doing de novo assembly")
m = denovo_assembly(preprocess_reads.out.full_len_reads, ref_genome)
m = denovo_assembly(full_len_reads, ref_genome)
} else {
build_minimap_index(ref_genome)
println("Doing reference based transcript analysis")
m = reference_assembly(build_minimap_index.out.index, ref_genome, preprocess_reads.out.full_len_reads)
m = reference_assembly(build_minimap_index.out.index, ref_genome, full_len_reads)
}
split_bam(m.bam)
@ -380,10 +394,10 @@ workflow pipeline {
software_versions,
workflow_params,
params.denovo,
pychopper_report,
summariseConcatReads.out.summary
.join(m.stats)
.join(run_gffcompare.out.gffcmp_dir)
.join(preprocess_reads.out.report)
.join(merge_gff_bundles.out.gff)
.toList().transpose().toList())
@ -395,19 +409,18 @@ workflow pipeline {
.join(seq_for_transcriptome_build))
if (use_ref_ann){
results = preprocess_reads.out.report
results = run_gffcompare.output.gffcmp_dir
.concat(
run_gffcompare.output.gffcmp_dir,
m.stats,
get_transcriptome.out.flatMap(map_sample_ids_cls))
m.stats,
get_transcriptome.out.flatMap(map_sample_ids_cls))
.map {it -> it[1]}
.concat(makeReport.out.report)
}
if (!use_ref_ann && !params.denovo){
results = preprocess_reads.out.report
.concat(m.stats,
get_transcriptome.out.flatMap(map_sample_ids_cls))
results = m.stats
.concat(
get_transcriptome.out.flatMap(map_sample_ids_cls))
.map {it -> it[1]}
.concat(makeReport.out.report)
@ -415,16 +428,17 @@ workflow pipeline {
if (params.denovo){
results = m.cds
.concat(m.stats,
seq_for_transcriptome_build,
get_transcriptome.out.flatMap(map_sample_ids_cls),
merge_gff_bundles.out.gff,
m.opt_qual_ch.flatMap {it ->
l = []
for (x in it[1..-1]){
l.add(tuple(it[0], x))
}
return l
})
seq_for_transcriptome_build,
get_transcriptome.out.flatMap(map_sample_ids_cls),
merge_gff_bundles.out.gff,
m.opt_qual_ch.flatMap {
it ->
l = []
for (x in it[1..-1]){
l.add(tuple(it[0], x))
}
return l
})
.map {it -> it[1]}
.concat(makeReport.out.report)
}

3
nextflow.config
View File

@ -33,7 +33,8 @@ params {
show_hidden_params = false
schema_ignore_params = 'show_hidden_params,validate_params,monochrome_logs,aws_queue,aws_image_prefix,wfversion'
// Process cDNA reads using pychopper, turn off for direct RNA:
direct_rna = true
// Options passed to pychopper:
pychopper_opts = "-m edlib"

5
nextflow_schema.json
View File

@ -66,6 +66,11 @@
"description": "Extra pychopper opts",
"default": "-m edlib"
},
"direct_rna": {
"type": "boolean",
"description": "Set to true for direct RNA sequencing. Omits the pychopper step.",
"default": false
},
"bundle_min_reads": {
"type": "integer",
"description": "Minimum size of bam bundle for parallel processing."