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

Resolve CW-3468 "Incorrect ids"

Browse Source
This commit is contained in:
Neil Horner 2024-02-13 08:48:21 +00:00
parent 4f671051c1
commit bb9acf55be
9 changed files with 79 additions and 128 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.pre-commit-config.yaml
View File

@ -8,7 +8,7 @@ repos:
always_run: true
pass_filenames: false
additional_dependencies:
- epi2melabs>=0.0.51
- epi2melabs>=0.0.52
- id: build_models
name: build_models
entry: datamodel-codegen --strict-nullable --base-class workflow_glue.results_schema_helpers.BaseModel --use-schema-description --disable-timestamp --input results_schema.yml --input-file-type openapi --output bin/workflow_glue/results_schema.py

2
CHANGELOG.md
View File

@ -4,7 +4,7 @@ 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]
## [v1.1.0]
### Changed
- Improve documentation around filtering of transcripts done before DTU analysis.
- Renamed files:

7
README.md
View File

@ -122,6 +122,13 @@ input_reads.fastq ─── input_directory ─── input_directory
| analyse_unclassified | boolean | Analyse unclassified reads from input directory. By default the workflow will not process reads in the unclassified directory. | If selected and if the input is a multiplex directory the workflow will also process the unclassified directory. | False |
### Output Options
| Nextflow parameter name | Type | Description | Help | Default |
|--------------------------|------|-------------|------|---------|
| out_dir | string | Directory for output of all user-facing files. | | output |
### Sample Options
| Nextflow parameter name | Type | Description | Help | Default |

109
bin/workflow_glue/report.py
View File

@ -51,23 +51,17 @@ def argparser():
"--gff_annotation", required=False, nargs='+',
help="transcriptome annotation gff file")
parser.add_argument(
"--gffcompare_dir", required=False, default=None, nargs='*',
"--gffcompare_dir", required=False, default=None,
help="gffcompare outout dir")
parser.add_argument(
"--pychop_report", required=False, default=None,
help="TSV summary file of pychopper statistics")
parser.add_argument(
"--sample_ids", required=True, nargs='+',
help="List of sample ids")
parser.add_argument(
"--isoform_table", required=False, type=Path,
help="Path to directory of TSV files with isoform summaries")
parser.add_argument(
"--isoform_table_nrows", required=False, type=int, default=5000,
help="Maximum rows to display in isoforms table")
parser.add_argument(
"--cluster_qc_dirs", required=False, type=str, default=None, nargs='*',
help="Directory with various cluster quality csvs")
parser.add_argument(
"--jaffal_csv", required=False, type=str, default=None,
help="Path to JAFFAL results csv")
@ -345,7 +339,7 @@ def grouped_bar(df, title="", tilted_xlabs=False):
return p
def gff_compare_plots(report, gffcompare_outdirs, sample_ids):
def gff_compare_plots(report, gffcompare_outdirs):
"""Create various sections and plots in a WfReport.
:param report: aplanat WFReport
@ -383,7 +377,10 @@ def gff_compare_plots(report, gffcompare_outdirs, sample_ids):
tabs = []
gff_fails = False
for id_, dir_ in zip(sample_ids, gffcompare_outdirs):
sample_ids = []
for dir_ in gffcompare_outdirs:
sample_id = dir_.name
sample_ids.append(sample_id) # Get sample ids fromt the folder name
stats, _, miss, novel, total = \
parse_gffcmp_stats(dir_ / 'str_merged.stats')
@ -396,7 +393,7 @@ def gff_compare_plots(report, gffcompare_outdirs, sample_ids):
tabs.append(Panel(
child=gridplot(
[bar_totals, bar_performance, bar_missed, bar_novel],
ncols=2, width=350, height=260), title=id_))
ncols=2, width=350, height=260), title=sample_id))
else:
gff_fails = True
@ -445,7 +442,7 @@ def gff_compare_plots(report, gffcompare_outdirs, sample_ids):
track_files = [x / 'str_merged.tracking' for x in gffcompare_outdirs]
df_tracking = load_sample_data(
df_tracking = load_data_add_sample_id(
track_files, sample_ids,
read_func=lambda x: pd.read_csv(
x, sep="\t", header=None,
@ -531,7 +528,7 @@ def gff_compare_plots(report, gffcompare_outdirs, sample_ids):
sys.stderr("Cannot find .tmap files in {}".format(gffcompare_outdirs))
return
df_tmap = load_sample_data(tmap_files, sample_ids)
df_tmap = load_data_add_sample_id(tmap_files, sample_ids)
for id_, df in df_tmap.groupby('sample_id'):
@ -596,59 +593,6 @@ def pychopper_plots(report, pychop_report):
section.plot(grid)
def cluster_quality(cluster_qc_dir, report, sample_ids):
"""Make cluster quality section."""
section = report.add_section()
section.markdown('''
### De novo clustering quality
This section shows plots relating to the clustering quality performed
by isONclust2. The full length reads are mapped to a reference genome
to create a ground truth of reads mapped to clusters. This is then compared
to the de novo-generated clusters, and the following statistics are
generated.
* [Homogeneity](https://scikit-learn.org/stable/modules/generated/
sklearn.metrics.homogeneity_score.html): Penalises over-clustering.
* [Completeness](https://scikit-learn.org/stable/modules/generated/
sklearn.metrics.completeness_score.html): Penalises under-clustering.
* [V-measure](https://clusteringjl.readthedocs.io/en/latest/vmeasure.html):
The harmonic mean of the homogeneity and completeness
* [Adjusted Rand Index](https://scikit-learn.org/stable/modules/generated/
sklearn.metrics.adjusted_rand_score.html): Intuitively, measures the
percentage of read pairs correctly clustered, normalized so that a perfect
clustering = 1 and a random cluster assignment achieves = 0
* NonSingleton: Clusters with multiple reads
* Singleton: Clusters consisting of a single read (These do not contribute
to the final transcript calling - I need to check this!)
''')
tabs = []
for id_, cluster_dir in zip(sample_ids, cluster_qc_dir):
plots = []
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)
bar = bars.simple_bar(
df.Statistic.values.tolist(), df.Value.values.tolist(),
colors=Colors.cerulean
)
bar.xaxis.major_label_orientation = math.pi / 2.8
plots.append(bar)
tabs.append(Panel(
child=gridplot(plots, ncols=4,
width=300, height=300), title=id_))
cover_panel = Tabs(tabs=tabs)
section.plot(cover_panel)
def transcript_table(report, isoform_table, max_rows):
"""Create searchable table of transcripts.
@ -685,7 +629,7 @@ def transcript_table(report, isoform_table, max_rows):
section.table(df, index=False)
def transcriptome_summary(report, gffs, sample_ids):
def transcriptome_summary(report, gffs):
"""
Plot transcriptome summaries.
@ -694,7 +638,6 @@ def transcriptome_summary(report, gffs, sample_ids):
:param report: aplanat WFReport
:param gffs: list of paths to gff transcriptome annotations
:param sample_ids: list of sample ids
"""
# test.db gets written to the git repo.
section = report.add_section()
@ -703,7 +646,8 @@ def transcriptome_summary(report, gffs, sample_ids):
''')
tabs = []
for id_, gff in zip(sample_ids, gffs):
for gff in gffs:
sample_id = Path(gff).name
plots = []
@ -747,7 +691,7 @@ def transcriptome_summary(report, gffs, sample_ids):
plots.append(bar_isos)
box = bars.boxplot_series(
[id_] * len(transcript_lens), transcript_lens,
[sample_id] * len(transcript_lens), transcript_lens,
width=70, ylim=(min(transcript_lens), max(transcript_lens)),
title='transcript lengths')
plots.append(box)
@ -780,14 +724,14 @@ def transcriptome_summary(report, gffs, sample_ids):
tabs.append(Panel(
child=gridplot(plots, ncols=4,
width=300, height=300), title=id_))
width=300, height=300), title=sample_id))
cover_panel = Tabs(tabs=tabs)
section.plot(cover_panel)
def load_sample_data(files, sample_ids, read_func=None):
"""Load CSVs into dataframe, and assign sample_id column."""
def load_data_add_sample_id(files, sample_ids, read_func=None):
"""Load CSVs and concat into single into dataframe, and assign sample_id column."""
df_ = pd.DataFrame()
if not files:
return None
@ -890,18 +834,18 @@ def de_section(report):
def main(args):
"""Run the entry point."""
sample_ids = args.sample_ids
sample_ids.sort()
report = WFReport(
"Transcript isoform report", "wf-transcriptomes",
revision=args.revision, commit=args.commit)
# QC
seq_stats_tabs(report, args.stats)
if args.alignment_stats is not None:
df_aln_stats = load_sample_data(args.alignment_stats, sample_ids)
stats_dfs = []
for stats_file in args.alignment_stats:
df = pd.read_csv(stats_file, sep='\t+')
stats_dfs.append(df)
aln_stats_df = pd.concat(stats_dfs)
section = report.add_section()
section.markdown('''
### Read mapping summary
@ -910,28 +854,23 @@ def main(args):
[seqkit](https://bioinf.shenwei.me/seqkit/)
`seqkit bam -s`''')
section.table(df_aln_stats)
section.table(aln_stats_df)
if args.pychop_report is not None:
pychopper_plots(report, args.pychop_report)
# Results
if args.gff_annotation is not None:
transcriptome_summary(
report, args.gff_annotation, sample_ids)
transcriptome_summary(report, args.gff_annotation)
if args.gffcompare_dir is not None:
gff_compare_plots(
report,
[Path(x) for x in args.gffcompare_dir],
sample_ids)
[x for x in Path(args.gffcompare_dir).iterdir()])
if args.isoform_table is not None:
transcript_table(report, args.isoform_table, args.isoform_table_nrows)
if args.cluster_qc_dirs is not None:
cluster_quality(args.cluster_qc_dirs, report, sample_ids)
if args.de_report:
de_section(report)

7
docs/06_input_parameters.md
View File

@ -11,6 +11,13 @@
| analyse_unclassified | boolean | Analyse unclassified reads from input directory. By default the workflow will not process reads in the unclassified directory. | If selected and if the input is a multiplex directory the workflow will also process the unclassified directory. | False |
### Output Options
| Nextflow parameter name | Type | Description | Help | Default |
|--------------------------|------|-------------|------|---------|
| out_dir | string | Directory for output of all user-facing files. | | output |
### Sample Options
| Nextflow parameter name | Type | Description | Help | Default |

71
main.nf
View File

@ -150,19 +150,18 @@ process preprocess_reads {
tuple val("${meta.alias}"),
path("${meta.alias}_pychopper_output/"),
emit: pychopper_output
path "${meta.alias}_pychopper_output/${meta.alias}_pychopper.tsv",
path("${meta.alias}_pychopper_output/pychopper.tsv"),
emit: report
script:
def cdna_kit = params.cdna_kit.split("-")[-1]
def extra_params = params.pychopper_opts ?: ''
def extra_params = params.pychopper_opts ?: ''
"""
pychopper -t ${params.threads} -k ${cdna_kit} -m ${params.pychopper_backend} ${extra_params} 'seqs.fastq.gz' ${meta.alias}_full_length_reads.fastq
mv pychopper.tsv ${meta.alias}_pychopper.tsv
workflow-glue generate_pychopper_stats --data ${meta.alias}_pychopper.tsv --output .
workflow-glue generate_pychopper_stats --data pychopper.tsv --output .
# Add sample id colum
sed "1s/\$/\tsample_id/; 1 ! s/\$/\t${meta.alias}/" ${meta.alias}_pychopper.tsv > tmp
mv tmp ${meta.alias}_pychopper.tsv
# Add sample id column
sed "1s/\$/\tsample_id/; 1 ! s/\$/\t${meta.alias}/" pychopper.tsv > tmp
mv tmp pychopper.tsv
mkdir "${meta.alias}_pychopper_output/"
@ -274,9 +273,9 @@ process merge_gff_bundles{
input:
tuple val(sample_id), path (gff_bundle)
output:
tuple val(sample_id), path('*.gff'), emit: gff
tuple val(sample_id), path("${sample_id}.gff"), emit: gff
script:
def merged_gff = "transcripts_${sample_id}.gff"
def merged_gff = "${sample_id}.gff"
"""
echo '##gff-version 2' >> $merged_gff;
echo '#pipeline-nanopore-isoforms: stringtie' >> $merged_gff;
@ -303,13 +302,12 @@ process run_gffcompare{
tuple val(sample_id), path(query_annotation)
path ref_annotation
output:
tuple val(sample_id), path("${sample_id}_gffcompare"),
emit: gffcmp_dir
tuple val(sample_id), path("${sample_id}"), emit: gffcmp_dir
path ("${sample_id}_annotated.gtf"), emit: gtf
tuple val(sample_id), path("${sample_id}_transcripts_table.tsv"),
emit: isoforms_table
script:
def out_dir = "${sample_id}_gffcompare"
def out_dir = "${sample_id}"
"""
mkdir $out_dir
echo "Doing comparison of reference annotation: ${ref_annotation} and the query annotation"
@ -320,14 +318,14 @@ process run_gffcompare{
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
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 \
--sample_id $sample_id \
--gffcompare_dir "${sample_id}_gffcompare"
--sample_id "${sample_id}" \
--gffcompare_dir "${out_dir}"
"""
}
@ -398,10 +396,9 @@ process makeReport {
path "params.json"
path "pychopper_report/*"
path"jaffal_csv/*"
val sample_ids
path "per_read_stats/?.gz"
path "aln_stats/*"
path gffcmp_dir
path "gffcmp_dir/*"
path "gff_annotation/*"
path "de_report/*"
path "seqkit/*"
@ -415,8 +412,6 @@ process makeReport {
path ("filtered_transcript_counts_with_genes.tsv"), emit: filtered, optional: true
path ("all_gene_counts.tsv"), emit: gene_counts, optional: true
shell:
// Convert the sample_id arrayList.
sids = new BlankSeparatedList(sample_ids)
report_name = "wf-transcriptomes-report.html"
'''
if [ -f "de_report/OPTIONAL_FILE" ]; then
@ -430,10 +425,10 @@ process makeReport {
else
OPT_GFF_ANNOTATION="--gff_annotation gff_annotation/*"
fi
if [ -f "OPTIONAL_FILE" ]; then
if [ -f "gffcmp_dir/OPTIONAL_FILE" ]; then
OPT_GFFCMP_DIR=""
else
OPT_GFFCMP_DIR="--gffcompare_dir !{gffcmp_dir}"
OPT_GFFCMP_DIR="--gffcompare_dir gffcmp_dir/"
fi
if [ -f "jaffal_csv/OPTIONAL_FILE" ]; then
OPT_JAFFAL_CSV=""
@ -460,7 +455,6 @@ process makeReport {
--params params.json \
${OPT_ALN} \
${OPT_PC_REPORT} \
--sample_ids !{sids} \
--stats per_read_stats/* \
${OPT_GFF_ANNOTATION} \
${OPT_ISO_TABLE} \
@ -596,7 +590,7 @@ workflow pipeline {
}
else{
full_len_reads = input_reads.map{ meta, reads -> [meta.alias, reads]}
pychopper_report = file("$projectDir/data/OPTIONAL_FILE")
pychopper_report = OPTIONAL_FILE
}
if (params.transcriptome_source != "precomputed"){
build_minimap_index(ref_genome)
@ -624,7 +618,7 @@ workflow pipeline {
optional_channel = Channel.fromPath("$projectDir/data/OPTIONAL_FILE")
gff_tuple = merge_gff_bundles.out.gff.combine(optional_channel)
gff_compare = OPTIONAL_FILE
isoforms_table = file("$projectDir/data/OPTIONAL_FILE")
isoforms_table = OPTIONAL_FILE
}
// For reference based assembly, there is only one reference
// So map this reference to all sample_ids
@ -638,10 +632,10 @@ workflow pipeline {
results = results.concat(assembly.bam.map {sample_id, bam, bai -> [bam, bai]}.flatten())
}
else{
gff_compare = file("$projectDir/data/OPTIONAL_FILE")
isoforms_table = file("$projectDir/data/OPTIONAL_FILE")
merge_gff = file("$projectDir/data/OPTIONAL_FILE")
assembly_stats = file("$projectDir/data/OPTIONAL_FILE")
gff_compare = OPTIONAL_FILE
isoforms_table = OPTIONAL_FILE
merge_gff = OPTIONAL_FILE
assembly_stats = OPTIONAL_FILE
use_ref_ann = false
}
@ -649,7 +643,7 @@ workflow pipeline {
gene_fusions(full_len_reads, jaffal_refBase, jaffal_genome, jaffal_annotation)
jaffal_out = gene_fusions.out.results_csv.map{ alias, csv -> csv}.collectFile(keepHeader: true, name: 'jaffal.csv')
} else{
jaffal_out = file("$projectDir/data/OPTIONAL_FILE")
jaffal_out = OPTIONAL_FILE
}
if (params.de_analysis){
@ -674,8 +668,8 @@ workflow pipeline {
de_outputs = de.de_outputs
counts = de.counts
} else{
de_report = file("$projectDir/data/OPTIONAL_FILE")
count_transcripts_file = file("$projectDir/data/OPTIONAL_FILE")
de_report = OPTIONAL_FILE
count_transcripts_file = OPTIONAL_FILE
}
makeReport(
@ -683,7 +677,6 @@ workflow pipeline {
workflow_params,
pychopper_report,
jaffal_out,
input_reads.map{ meta, fastq -> meta.alias}.collect(),
per_read_stats,
assembly_stats,
gff_compare,
@ -724,8 +717,8 @@ workflow pipeline {
if (params.de_analysis){
de_results = report.concat(
transcriptome, de_outputs.flatten(), counts.flatten(),
makeReport.out.results_dge, makeReport.out.tpm,
makeReport.out.filtered, makeReport.out.unfiltered,
makeReport.out.results_dge, makeReport.out.tpm,
makeReport.out.filtered, makeReport.out.unfiltered,
makeReport.out.gene_counts)
// Output de_analysis results in the dedicated directory.
results = results.concat(de_results.map{ [it, "de_analysis"] })
@ -768,7 +761,7 @@ workflow {
error = "--ref_genome: File doesn't exist, check path."
}
}else {
ref_genome = file("$projectDir/data/OPTIONAL_FILE")
ref_genome = OPTIONAL_FILE
}
if (params.containsValue("denovo")) {
error = "Denovo transcriptome source is no longer supported. Please use the reference-guided or precomputed options."
@ -782,7 +775,7 @@ workflow {
}
use_ref_ann = true
}else{
ref_annotation= file("$projectDir/data/OPTIONAL_FILE")
ref_annotation= OPTIONAL_FILE
use_ref_ann = false
}
if (params.jaffal_refBase){
@ -793,7 +786,7 @@ workflow {
}else{
jaffal_refBase = null
}
ref_transcriptome = file("$projectDir/data/OPTIONAL_FILE")
ref_transcriptome = OPTIONAL_FILE
if (params.ref_transcriptome){
log.info("Reference Transcriptome provided will be used for differential expression.")
ref_transcriptome = file(params.ref_transcriptome, type:"file")

2
nextflow.config
View File

@ -116,7 +116,7 @@ manifest {
description = 'Transcriptome analysis including gene fusions, differential expression as well as assembly and annotation of cDNA and direct RNA sequencing data.'
mainScript = 'main.nf'
nextflowVersion = '>=23.04.2'
version = 'v1.0.0'
version = 'v1.1.0'
}
epi2melabs {

2
subworkflows/JAFFAL/gene_fusions.nf
View File

@ -17,7 +17,7 @@ process jaffal{
# JAFFAL exists with status code 1 when there's 0 fusion hits. Prevent this with '||:'
$params.jaffal_dir/tools/bin/bpipe run \
-n $params.threads \
-n "${task.cpus}" \
-p jaffa_output="\$JAFFAOUT/" \
-p refBase=$refBase \
-p genome=$genome \

5
subworkflows/reference_assembly.nf
View File

@ -30,6 +30,11 @@ process map_reads{
| 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
# Add sample id header and column
sed "s/\$/${sample_id}/" "${sample_id}_read_aln_stats.tsv" \
| sed "1 s/${sample_id}/sample_id/" > tmp
mv tmp "${sample_id}_read_aln_stats.tsv"
if [[ -s "internal_priming_fail.tsv" ]];
then
tail -n +2 "internal_priming_fail.tsv" | awk '{print ">" \$1 "\\n" \$4 }' - > "context_internal_priming_fail_start.fasta"