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

Merge branch 'ref-prep' into 'dev'

Browse Source 
[CW-7166] [CW-7169] Reference and annotation preparation

See merge request epi2melabs/workflows/wf-transcriptomes!239
This commit is contained in:
Kiah McIntosh 2026-05-07 10:09:36 +00:00
commit 9ddec6ff61
11 changed files with 1132 additions and 304 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

178
.gitlab-ci.yml
View File

@ -56,15 +56,9 @@ docker-run:
parallel:
matrix:
- MATRIX_NAME: [
"differential_expression", "isoforms", "isoforms_bam",
"only_differential_expression", "differential_expression_gff3",
"ncbi_gzip", "ncbi_no_gene_id", "ensembl_with_versions",
"differential_expression_mouse", "no_ref_annotation",
"unstranded_annotation_warn", "igv", "igv_fai_gz",
"mismatch_seqid_error","mismatch_seqid_warn",
"codex_discover", "codex_fixed_annotation", "codex_direct_rna",
"codex_de_analysis", "codex_invalid_mode",
"codex_conflicting_library_flags"
"de", "discover", "igv",
"smoke_discover", "smoke_fixed", "smoke_direct_rna", "smoke_de",
"no_annotation", "invalid_mode", "conflicting_flags"
]
rules:
# NOTE As we're overriding the rules block for the included docker-run
@ -72,20 +66,13 @@ docker-run:
# being incorrectly scheduled for "detached merge request pipelines" etc.
- if: ($CI_COMMIT_BRANCH == null || $CI_COMMIT_BRANCH == "dev-template")
when: never
- if: $MATRIX_NAME == "isoforms"
- if: $MATRIX_NAME == "discover"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq \
--ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa --ref_annotation ${CI_PROJECT_NAME}/data/chr20/gencode.v22.annotation.chr20.gtf --cdna_preprocess --pychopper_backend phmm"
NF_IGNORE_PROCESSES: preprocess_reads,filter_unstranded_annotation,validate_ref_annotation,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
- if: $MATRIX_NAME == "isoforms_bam"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS: "--bam ${CI_PROJECT_NAME}/data/ERR6053095_chr20.bam \
--ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa --ref_annotation ${CI_PROJECT_NAME}/data/chr20/gencode.v22.annotation.chr20.gtf --cdna_preprocess --pychopper_backend phmm \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: preprocess_reads,filter_unstranded_annotation,validate_ref_annotation,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
- if: $MATRIX_NAME == "no_ref_annotation"
- if: $MATRIX_NAME == "no_annotation"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq \
@ -94,7 +81,7 @@ docker-run:
NF_IGNORE_PROCESSES: preprocess_reads
ASSERT_NEXTFLOW_FAILURE: "1"
ASSERT_NEXTFLOW_FAILURE_REXP: "Missing required parameter: --ref_annotation"
- if: $MATRIX_NAME == "differential_expression"
- if: $MATRIX_NAME == "de"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
@ -105,101 +92,6 @@ docker-run:
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: preprocess_reads,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
- if: $MATRIX_NAME == "only_differential_expression"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
--transcriptome_mode fixed_annotation \
--de_analysis \
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression/hg38_chr20.fa \
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.gff \
--direct_rna --minimap2_index_opts '-k 15' \
--sample_sheet test_data/sample_sheet.csv \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: >
preprocess_reads,faidx,validate_ref_annotation,gz_faidx,merge_transcriptomes,merge_gff_bundles,assemble_transcripts,decompress_annotation,decompress_ref,
build_minimap_index,get_transcriptome,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
- if: $MATRIX_NAME == "differential_expression_gff3"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
--transcriptome_mode fixed_annotation \
--de_analysis \
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression/hg38_chr20.fa \
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.gff3 \
--direct_rna --minimap2_index_opts '-k 15' \
--sample_sheet test_data/sample_sheet.csv \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: >
preprocess_reads,validate_ref_annotation,faidx,gz_faidx,merge_transcriptomes,assemble_transcripts,decompress_annotation,decompress_ref,
build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
- if: $MATRIX_NAME == "ncbi_gzip"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression_ncbi.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS:
"--fastq ${CI_PROJECT_NAME}/data/differential_expression_ncbi/differential_expression_fastq \
--de_analysis \
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GRCh38.p14.NCBI_test.fna.gz \
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GRCh38.p14_NCBI_test.gtf.gz \
--direct_rna --minimap2_index_opts '-w 25' \
--sample_sheet test_data/sample_sheet.csv \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: >
preprocess_reads,faidx,gz_faidx,merge_transcriptomes,assemble_transcripts,
build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
- if: $MATRIX_NAME == "ncbi_no_gene_id"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression_ncbi.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS:
"--fastq ${CI_PROJECT_NAME}/data/differential_expression_ncbi/differential_expression_fastq \
--transcriptome_mode fixed_annotation --de_analysis \
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.fna.gz \
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.gff.gz \
--direct_rna \
--sample_sheet test_data/sample_sheet.csv \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: >
preprocess_reads,faidx,validate_ref_annotation,gz_faidx,merge_transcriptomes,assemble_transcripts,
build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
- if: $MATRIX_NAME == "ensembl_with_versions"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS:
"--fastq ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
--de_analysis \
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz \
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/Homo_sapiens.GRCh38.109.gtf.gz \
--direct_rna --transcriptome_mode fixed_annotation \
--sample_sheet test_data/sample_sheet.csv \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: >
preprocess_reads,validate_ref_annotation,faidx,gz_faidx,merge_transcriptomes,assemble_transcripts,
build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
- if: $MATRIX_NAME == "differential_expression_mouse"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression_mouse.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression_mouse.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression_mouse.tar.gz -C ${CI_PROJECT_NAME}/data/&& wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS:
"--fastq ${CI_PROJECT_NAME}/data/differential_expression_mouse/differential_expression_fastq \
--transcriptome_mode fixed_annotation --de_analysis \
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression_mouse/GRCm39.genome.fa.gz \
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression_mouse/gencode.vM33.annotation.gtf \
--direct_rna \
--sample_sheet ${CI_PROJECT_NAME}/data/differential_expression_mouse/sample_sheet.csv \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: >
preprocess_reads,validate_ref_annotation,faidx,gz_faidx,merge_transcriptomes,assemble_transcripts,decompress_annotation,
build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam
- if: $MATRIX_NAME == "unstranded_annotation_warn"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
--de_analysis \
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression/hg38_chr20.fa \
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/unstranded_annotation.gtf \
--direct_rna --minimap2_index_opts '-k 15' --sample_sheet ${CI_PROJECT_NAME}/data/differential_expression/sample_sheet.csv \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: preprocess_reads,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
AFTER_NEXTFLOW_CMD: "grep 'Unstranded entries found and excluded from differential expression analysis.' '.nextflow.log'"
- if: $MATRIX_NAME == "igv"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS: "--fastq ${CI_PROJECT_NAME}/data/differential_expression/differential_expression_fastq \
@ -213,53 +105,7 @@ docker-run:
NF_IGNORE_PROCESSES: >
preprocess_reads,gz_faidx,merge_transcriptomes,assemble_transcripts,decompress_annotation,decompress_ref,
build_minimap_index,validate_ref_annotation,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
- if: $MATRIX_NAME == "igv_fai_gz"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/differential_expression_ncbi.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS:
"--fastq ${CI_PROJECT_NAME}/data/differential_expression_ncbi/differential_expression_fastq \
--transcriptome_mode fixed_annotation --de_analysis \
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.fna.gz \
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.gff.gz \
--direct_rna \
--sample_sheet test_data/sample_sheet.csv \
--igv \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: >
preprocess_reads,validate_ref_annotation,merge_transcriptomes,assemble_transcripts,
build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome,faidx,gz_faidx
- if: $MATRIX_NAME == "mismatch_seqid_error"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression_ncbi.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS:
"--fastq ${CI_PROJECT_NAME}/data/differential_expression_ncbi/differential_expression_fastq \
--de_analysis \
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/hg38_chr20.fa \
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GRCh38.p14_NCBI_test.gtf.gz \
--direct_rna \
--sample_sheet test_data/sample_sheet.csv \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: >
preprocess_reads,faidx,gz_faidx,merge_transcriptomes,assemble_transcripts,
build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
ASSERT_NEXTFLOW_FAILURE: 1
ASSERT_NEXTFLOW_FAILURE_REXP: No overlapping seqnames were found between the reference annotation and the reference genome
- if: $MATRIX_NAME == "mismatch_seqid_warn"
variables:
NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression_ncbi.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression_ncbi.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -nv -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config;"
NF_WORKFLOW_OPTS:
"--fastq ${CI_PROJECT_NAME}/data/differential_expression_ncbi/differential_expression_fastq \
--de_analysis \
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GRCh38.p14.NCBI_test.fna.gz \
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GRCh38.p14_NCBI_test.gtf.gz \
--direct_rna \
--sample_sheet test_data/sample_sheet.csv --minimap2_index_opts '-w 25' \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config "
NF_IGNORE_PROCESSES: >
preprocess_reads,faidx,gz_faidx,merge_transcriptomes,assemble_transcripts,
build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam,decompress_transcriptome,preprocess_ref_transcriptome
AFTER_NEXTFLOW_CMD: "grep 'Some seqnames are present in the genome but not the annotation' '.nextflow.log'"
- if: $MATRIX_NAME == "codex_discover"
- if: $MATRIX_NAME == "smoke_discover"
variables:
NF_BEFORE_SCRIPT: ":"
NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf"
@ -268,7 +114,7 @@ docker-run:
test -f ${CI_PROJECT_NAME}/cohort/transcripts.gtf &&
test -f ${CI_PROJECT_NAME}/cohort/cohort.transcriptome.fa &&
test -f ${CI_PROJECT_NAME}/samples/sampleA/transcripts.gtf
- if: $MATRIX_NAME == "codex_fixed_annotation"
- if: $MATRIX_NAME == "smoke_fixed"
variables:
NF_BEFORE_SCRIPT: ":"
NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --transcriptome_mode fixed_annotation"
@ -276,7 +122,7 @@ docker-run:
AFTER_NEXTFLOW_CMD: >
test -f ${CI_PROJECT_NAME}/cohort/transcripts.gtf &&
test -f ${CI_PROJECT_NAME}/cohort/transcript_counts.tsv
- if: $MATRIX_NAME == "codex_direct_rna"
- if: $MATRIX_NAME == "smoke_direct_rna"
variables:
NF_BEFORE_SCRIPT: ":"
NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --direct_rna"
@ -284,7 +130,7 @@ docker-run:
AFTER_NEXTFLOW_CMD: >
test -f ${CI_PROJECT_NAME}/cohort/alignments/sampleA.aligned.sorted.bam &&
test -f ${CI_PROJECT_NAME}/cohort/sqanti_cohort/classification_summary.tsv
- if: $MATRIX_NAME == "codex_de_analysis"
- if: $MATRIX_NAME == "smoke_de"
variables:
NF_BEFORE_SCRIPT: ":"
NF_WORKFLOW_OPTS: "--fastq test_data/smoke/de --sample_sheet test_data/smoke/sample_sheet_de.csv --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --de_analysis --reference_level control --covariates batch"
@ -292,14 +138,14 @@ docker-run:
AFTER_NEXTFLOW_CMD: >
test -f ${CI_PROJECT_NAME}/de_analysis/condition_treated_vs_control/results_dge.tsv &&
test -f ${CI_PROJECT_NAME}/de_analysis/condition_treated_vs_control/results_dtu_transcript.tsv
- if: $MATRIX_NAME == "codex_invalid_mode"
- if: $MATRIX_NAME == "invalid_mode"
variables:
NF_BEFORE_SCRIPT: ":"
NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --transcriptome_mode nonsense"
NF_IGNORE_PROCESSES: preprocess_reads
ASSERT_NEXTFLOW_FAILURE: "1"
ASSERT_NEXTFLOW_FAILURE_REXP: "nonsense is not a valid choice"
- if: $MATRIX_NAME == "codex_conflicting_library_flags"
- if: $MATRIX_NAME == "conflicting_flags"
variables:
NF_BEFORE_SCRIPT: ":"
NF_WORKFLOW_OPTS: "--fastq test_data/smoke/reads.fastq --sample sampleA --ref_genome test_data/smoke/reference.fa --ref_annotation test_data/smoke/annotation.gtf --direct_rna --cdna_preprocess"

2
README.md
View File

@ -377,6 +377,8 @@ Output files may be aggregated including information for all samples or provided
| Aligned BAM | cohort/alignments/{{ alias }}.aligned.sorted.bam | Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV. | per-sample |
| Aligned BAM index | cohort/alignments/{{ alias }}.aligned.sorted.bam.bai | Index for the aligned BAM. | per-sample |
| Alignment summary | cohort/alignments/{{ alias }}.flagstat.txt | samtools flagstat output for the aligned BAM. | per-sample |
| Reference and annotation preparation summary | cohort/reference/annotation_reference_summary.json | Summary of reference and annotation preparation, including seqname overlap, build/provider hints, and excluded unstranded annotation counts. | aggregated |
| Excluded unstranded annotation records | cohort/reference/unstranded_annotation.gtf | Full set of annotation records excluded because their strand was not '+' or '-'. Present only when unstranded records are found. | aggregated |
| Cohort transcriptome GTF | cohort/transcripts.gtf | Joint bambu transcript model used as the primary cohort transcriptome. | aggregated |
| Cohort transcriptome FASTA | cohort/cohort.transcriptome.fa | Transcript sequences derived from the joint cohort GTF. | aggregated |
| Cohort transcript counts | cohort/transcript_counts.tsv | Transcript-level count matrix produced by bambu. | aggregated |

533
bin/workflow_glue/prepare_annotation_reference.py Normal file
View File

@ -0,0 +1,533 @@
"""Prepare reference and annotation files for transcriptome analysis."""
import gzip
import json
from pathlib import Path
import re
import shutil
import subprocess
import sys
from .util import get_named_logger, wf_parser # noqa: ABS101
UNSTRANDED_WARNING = (
"Warning: Unstranded entries found and excluded from downstream "
"transcriptome analysis."
)
BUILD_PATTERNS = {
"GRCh38": re.compile(r"(?<![A-Za-z0-9])(?:GRCh38|hg38)(?![A-Za-z0-9])", re.I),
"GRCh37": re.compile(r"(?<![A-Za-z0-9])(?:GRCh37|hg19)(?![A-Za-z0-9])", re.I),
"GRCm39": re.compile(r"(?<![A-Za-z0-9])GRCm39(?![A-Za-z0-9])", re.I),
"GRCm38": re.compile(r"(?<![A-Za-z0-9])GRCm38(?![A-Za-z0-9])", re.I),
}
PROVIDER_PATTERNS = {
"GENCODE": re.compile(r"(?<![A-Za-z0-9])gencode(?![A-Za-z0-9])", re.I),
"Ensembl": re.compile(r"(?<![A-Za-z0-9])ensembl(?![A-Za-z0-9])", re.I),
"RefSeq": re.compile(r"(?<![A-Za-z0-9])(?:refseq|ncbi)(?![A-Za-z0-9])", re.I),
}
def _is_gzip(path):
"""Return whether path looks gzip-compressed."""
return Path(path).suffix.lower() == ".gz"
def _is_gff(path):
"""Return whether path looks like a GFF/GFF3 annotation."""
p = Path(path)
# check all suffixes for GFF extensions
suffixes = {s.lower() for s in p.suffixes}
return bool(suffixes & {".gff", ".gff3"})
def _find_patterns(text, patterns):
"""Find matching pattern names in text."""
return {name for name, pattern in patterns.items() if pattern.search(text)}
class PreparedReference:
"""A reference genome prepared for analysis."""
def __init__(self, input_path, work_dir):
"""Initialize and prepare reference genome."""
self.input_path = Path(input_path)
self.work_dir = Path(work_dir)
self.output_path = self.work_dir / "reference.fasta"
self._seqnames = None
self._prepare()
def _prepare(self):
"""Prepare reference (decompress if gzipped)."""
if _is_gzip(self.input_path):
with gzip.open(self.input_path, "rb") as src, open(
self.output_path, "wb"
) as dst:
shutil.copyfileobj(src, dst)
else:
shutil.copyfile(self.input_path, self.output_path)
if not self.output_path.exists() or self.output_path.stat().st_size < 1:
raise ValueError(f"Prepared file is empty: {self.output_path}")
@property
def seqnames(self):
"""Extract and cache sorted FASTA sequence names."""
if self._seqnames is None:
ids = set()
with open(self.output_path, encoding="utf-8") as f:
for line in f:
if line.startswith(">"):
ids.add(line[1:].split()[0])
self._seqnames = sorted(ids)
return self._seqnames
def detect_hints(self):
"""Detect build and provider hints from filename and headers."""
builds = set()
providers = set()
# check filename
builds.update(_find_patterns(self.input_path.name, BUILD_PATTERNS))
providers.update(_find_patterns(self.input_path.name, PROVIDER_PATTERNS))
# check headers
with open(self.output_path, encoding="utf-8") as f:
for line in f:
if line.startswith(">"):
builds.update(_find_patterns(line, BUILD_PATTERNS))
providers.update(_find_patterns(line, PROVIDER_PATTERNS))
return {
"builds": sorted(builds),
"providers": sorted(providers),
}
class AttributeSanitizer:
"""Fixes common GTF attribute problems."""
ATTRIBUTE_RE = re.compile(r'(^|;\s*){key}\s+"([^"]*)";?')
@classmethod
def _pattern_for(cls, key):
"""Return regex pattern for a GTF attribute key."""
return re.compile(cls.ATTRIBUTE_RE.pattern.format(key=re.escape(key)))
@classmethod
def _extract(cls, attr_field, key):
"""Extract a quoted GTF attribute value."""
match = cls._pattern_for(key).search(attr_field)
if not match:
return None
value = match.group(2)
# normalize whitespace and remove quotes/semicolons
value = re.sub(r"\s+", " ", value.replace('"', "").replace(";", "")).strip()
return value or None
@classmethod
def _replace(cls, attr_field, key, value):
"""Replace a quoted GTF attribute value."""
pattern = cls._pattern_for(key)
return pattern.sub(
lambda m: f'{m.group(1)}{key} "{value}"; ', # noqa: E702
attr_field,
count=1,
)
@classmethod
def fix(cls, attr_field):
"""Fix gene_id if missing or set to literal 'transcript_id'."""
transcript_id = cls._extract(attr_field, "transcript_id")
if transcript_id is None:
return attr_field, False
gene_id = cls._extract(attr_field, "gene_id")
# add gene_id if missing
if gene_id is None:
return (
f'gene_id "{transcript_id}"; {attr_field.strip()}', # noqa: E702
True,
)
# fix gene_id="transcript_id" literal
if gene_id == "transcript_id":
return cls._replace(attr_field, "gene_id", transcript_id), True
return attr_field, False
class Annotation:
"""An annotation file with preparation and filtering capabilities."""
MAX_UNSTRANDED_EXAMPLES = 20
def __init__(self, input_path, work_dir, gffread="gffread"):
"""Initialize and convert annotation to GTF."""
self.input_path = Path(input_path)
self.work_dir = Path(work_dir)
self.is_gff = _is_gff(input_path)
self.gffread = gffread
# outputs
self.unfiltered_path = None
self.prepared_path = None
self._convert_to_gtf()
def _open_gtf(self, path):
"""Open a GTF file, handling gzip transparently."""
if _is_gzip(path):
return gzip.open(path, "rt", encoding="utf-8")
return open(path, encoding="utf-8")
def _read_headers(self, path):
"""Yield comment/header lines from a GTF file."""
with self._open_gtf(path) as f:
for line in f:
if line.startswith("#"):
yield line.rstrip("\n")
def _read_records(self, path):
"""Yield (line_text, fields_list) for all lines in a GTF file."""
with self._open_gtf(path) as f:
for line in f:
stripped = line.rstrip("\n")
if not stripped or stripped.startswith("#"):
yield line, None
else:
fields = stripped.split("\t")
if len(fields) != 9:
raise ValueError(
f"Malformed annotation record has "
f"{len(fields)} columns; expected 9 GTF " # noqa: E702
f"columns: {stripped}"
)
yield line, fields
def _extract_seqnames(self, path):
"""Extract sorted unique seqnames from a GTF file."""
ids = set()
with self._open_gtf(path) as f:
for line in f:
stripped = line.rstrip("\n")
if stripped and not stripped.startswith("#"):
fields = stripped.split("\t")
if len(fields) >= 1:
ids.add(fields[0])
return sorted(ids)
def _convert_to_gtf(self):
"""Convert GFF/GFF3 to GTF using gffread, or use input if already GTF."""
if self.is_gff:
dest = self.work_dir / "annotation_unfiltered.gtf"
if _is_gzip(self.input_path):
# decompress as gff doesn't support gzip input
name_without_gz = self.input_path.name[:-3]
suffix = Path(name_without_gz).suffix or ".gff"
intermediate = self.work_dir / f"annotation_input{suffix}"
with gzip.open(self.input_path, "rb") as src, open(
intermediate, "wb"
) as dst:
shutil.copyfileobj(src, dst)
else:
intermediate = self.input_path
# use gffread to convert to GTF
result = subprocess.run(
[self.gffread, "-T", str(intermediate), "-o", str(dest)],
check=False,
capture_output=True,
text=True,
)
if result.returncode != 0:
details = "\n".join(
part
for part in (result.stdout, result.stderr)
if part
)
raise ValueError(
"Failed to convert annotation to GTF with "
"gffread:\n" + details
)
if not dest.exists() or dest.stat().st_size < 1:
raise ValueError(f"Prepared annotation is empty: {dest}")
self.unfiltered_path = dest
else:
# already GTF
if not self.input_path.exists() or self.input_path.stat().st_size < 1:
raise ValueError(f"Prepared annotation is empty: {self.input_path}")
self.unfiltered_path = self.input_path
@property
def seqnames(self):
"""Get seqnames from prepared (stranded) annotation."""
if self.prepared_path is None:
raise ValueError("Must call filter_stranded() before accessing seqnames")
return self._extract_seqnames(self.prepared_path)
def filter_stranded(self):
"""Filter to stranded records, return statistics."""
dest = self.work_dir / "annotation.gtf"
unstranded_dest = self.work_dir / "unstranded_annotation.gtf"
stats = {
"total_records": 0,
"kept_records": 0,
"excluded_unstranded_records": 0,
"sanitised_attribute_records": 0,
"unstranded_examples": [],
}
unstranded_file = None
try:
with open(dest, "w", encoding="utf-8") as out:
for line, fields in self._read_records(self.unfiltered_path):
# pass through headers and blank lines
if fields is None:
out.write(line)
continue
stats["total_records"] += 1
if fields[6] in {"+", "-"}:
# fix attributes if needed of stranded records
fixed_attrs, was_fixed = AttributeSanitizer.fix(fields[8])
if was_fixed:
stats["sanitised_attribute_records"] += 1
fields[8] = fixed_attrs
out.write("\t".join(fields) + "\n")
stats["kept_records"] += 1
else:
# collect unstranded records for posterity
stats["excluded_unstranded_records"] += 1
if (
len(stats["unstranded_examples"])
< self.MAX_UNSTRANDED_EXAMPLES
):
stats["unstranded_examples"].append(
line.rstrip("\n")
)
if unstranded_file is None:
unstranded_file = open(
unstranded_dest, "w", encoding="utf-8"
)
unstranded_file.write(line)
finally:
if unstranded_file:
unstranded_file.close()
# clean up empty unstranded file
if stats["excluded_unstranded_records"] == 0 and unstranded_dest.exists():
unstranded_dest.unlink()
self.prepared_path = dest
stats["unstranded_path"] = (
str(unstranded_dest) if stats["excluded_unstranded_records"] > 0 else None
)
return stats
def detect_hints(self):
"""Detect build and provider hints from filename and headers."""
builds = set()
providers = set()
# check filename
builds.update(_find_patterns(self.input_path.name, BUILD_PATTERNS))
providers.update(_find_patterns(self.input_path.name, PROVIDER_PATTERNS))
# check headers
for line in self._read_headers(self.input_path):
builds.update(_find_patterns(line, BUILD_PATTERNS))
providers.update(_find_patterns(line, PROVIDER_PATTERNS))
return {
"builds": sorted(builds),
"providers": sorted(providers),
}
def validate_seqname_overlap(annotation_seqnames, reference_seqnames):
"""Return overlap statistics for annotation/reference seqnames."""
ann_set = set(annotation_seqnames)
ref_set = set(reference_seqnames)
matches = sorted(ann_set & ref_set)
only_in_annotation = sorted(ann_set - ref_set)
only_in_reference = sorted(ref_set - ann_set)
return {
"has_overlap": bool(matches),
"matches": matches,
"only_in_annotation": only_in_annotation,
"only_in_reference": only_in_reference,
}
def build_warnings(filter_stats, seqnames, ref_hints, ann_hints):
"""Build all warning messages."""
warnings = []
if filter_stats["excluded_unstranded_records"]:
warnings.extend([
UNSTRANDED_WARNING,
(
f"Excluded {filter_stats['excluded_unstranded_records']} "
f"unstranded annotation records; " # noqa: E702
f"kept {filter_stats['kept_records']} stranded records."
),
"A sample of unstranded entries:",
*filter_stats["unstranded_examples"],
])
if filter_stats["sanitised_attribute_records"]:
warnings.append(
f"Warning: Sanitised {filter_stats['sanitised_attribute_records']} "
"annotation records with missing or malformed gene_id attributes."
)
if seqnames["only_in_annotation"]:
warnings.extend([
"Warning: Some seqnames are present in the annotation but not the genome:",
*seqnames["only_in_annotation"][:5],
])
if seqnames["only_in_reference"]:
warnings.extend([
"Warning: Some seqnames are present in the genome but not the annotation:",
*seqnames["only_in_reference"][:5],
])
ref_builds = ref_hints["builds"]
ann_builds = ann_hints["builds"]
if (
ref_builds
and ann_builds
and not (set(ref_builds) & set(ann_builds))
):
warnings.append(
"Warning: Reference and annotation contain different genome "
f"build hints: reference={', '.join(ref_builds)}; " # noqa: E702
f"annotation={', '.join(ann_builds)}."
)
ref_providers = ref_hints["providers"]
ann_providers = ann_hints["providers"]
if (
ref_providers
and ann_providers
and not (set(ref_providers) & set(ann_providers))
):
warnings.append(
"Warning: Reference and annotation contain different "
f"provider hints: " # noqa: E702
f"reference={', '.join(ref_providers)}; " # noqa: E702
f"annotation={', '.join(ann_providers)}."
)
return warnings
def prepare_annotation_reference(annotation, reference, out_dir, gffread="gffread"):
"""Prepare annotation.gtf and reference.fasta in out_dir."""
out_dir = Path(out_dir)
out_dir.mkdir(parents=True, exist_ok=False)
# prepare reference genome
ref = PreparedReference(reference, out_dir)
# prepare annotation (convert to GTF if needed, then filter to stranded)
ann = Annotation(annotation, out_dir, gffread=gffread)
filter_stats = ann.filter_stranded()
# check seqname overlap between prepared annotation and reference
seqnames = validate_seqname_overlap(ann.seqnames, ref.seqnames)
if not seqnames["has_overlap"]:
raise ValueError(
"ERROR: No overlapping seqnames were found between the "
"reference annotation and the reference genome.\n"
"Annotation ID examples:\n"
+ "\n".join(seqnames["only_in_annotation"][:5])
+ "\n"
+ "Reference ID examples:\n"
+ "\n".join(seqnames["only_in_reference"][:5])
)
# heuristically detect build and provider hints from filenames and headers
ref_hints = ref.detect_hints()
ann_hints = ann.detect_hints()
warnings = build_warnings(filter_stats, seqnames, ref_hints, ann_hints)
# build a summary
summary = {
"annotation": {
"input": str(ann.input_path),
"prepared": str(ann.prepared_path),
"unfiltered": str(ann.unfiltered_path),
"was_gff": ann.is_gff,
**filter_stats,
},
"reference": {
"input": str(ref.input_path),
"prepared": str(ref.output_path),
},
"seqnames": seqnames,
"analysis_seqnames": seqnames,
"seqname_overlap": seqnames["matches"],
"only_in_annotation": seqnames["only_in_annotation"],
"only_in_reference": seqnames["only_in_reference"],
"reference_build_hints": ref_hints["builds"],
"annotation_build_hints": ann_hints["builds"],
"reference_provider_hints": ref_hints["providers"],
"annotation_provider_hints": ann_hints["providers"],
"warnings": warnings,
}
summary_path = out_dir / "annotation_reference_summary.json"
with open(summary_path, "w", encoding="utf-8") as f:
json.dump(summary, f, indent=2)
f.write("\n")
return summary
def main(args):
"""Run the annotation/reference preparation command."""
logger = get_named_logger("prepRef")
summary = prepare_annotation_reference(
args.annotation,
args.reference,
args.out_dir,
gffread=args.gffread,
)
for warning in summary["warnings"]:
sys.stdout.write(f"{warning}\n")
logger.info(
"Prepared annotation %s and reference %s.",
summary["annotation"]["prepared"],
summary["reference"]["prepared"],
)
def argparser():
"""Argument parser for the preparation entry point."""
parser = wf_parser("prepare_annotation_reference")
parser.add_argument(
"--annotation", required=True, help="Reference annotation GTF/GFF/GFF3."
)
parser.add_argument("--reference", required=True, help="Reference genome FASTA.")
parser.add_argument(
"--out_dir",
required=True,
help="Output directory for prepared annotation/reference files.",
)
parser.add_argument("--gffread", default="gffread", help="gffread executable.")
return parser

103
bin/workflow_glue/report.py
View File

@ -122,6 +122,22 @@ def _load_de_qc(de_dir):
return None
def _load_annotation_reference_summary(cohort_dir):
"""Load reference/annotation preparation summary JSON."""
summary_file = Path(cohort_dir) / "reference" / "annotation_reference_summary.json"
if summary_file.exists():
with open(summary_file) as f:
return json.load(f)
return None
def _format_hint_values(hints):
"""Format provenance hints for a compact table cell."""
if not hints:
return "None detected"
return ", ".join(hints)
def _create_warning_banner(message, level="warning"):
"""Create a styled warning banner."""
colors = {
@ -200,6 +216,93 @@ def main(args):
# Load bambu QC statistics
bambu_qc = _load_bambu_qc(args.cohort_dir)
annotation_reference_summary = _load_annotation_reference_summary(args.cohort_dir)
if annotation_reference_summary:
with report.add_section("Reference and Annotation Checks", "Reference"):
for warning in annotation_reference_summary.get("warnings", []):
_create_warning_banner(warning, level="warning")
seqname_rows = [
(
"Overlapping seqnames",
len(annotation_reference_summary.get("seqname_overlap", [])),
),
(
"Seqnames only in annotation",
len(annotation_reference_summary.get("only_in_annotation", [])),
),
(
"Seqnames only in reference",
len(annotation_reference_summary.get("only_in_reference", [])),
),
]
annotation_summary = annotation_reference_summary.get("annotation", {})
seqname_rows.extend([
(
"Annotation records retained",
annotation_summary.get("kept_records", "N/A"),
),
(
"Unstranded records excluded",
annotation_summary.get("excluded_unstranded_records", "N/A"),
),
(
"Annotation attributes sanitised",
annotation_summary.get("sanitised_attribute_records", "N/A"),
),
])
DataTable.from_pandas(
pd.DataFrame(seqname_rows, columns=["Check", "Value"]),
paging=False,
searchable=False,
)
with h3("Build and Provider Hints"):
hint_rows = [
(
"Reference build",
_format_hint_values(
annotation_reference_summary.get(
"reference_build_hints", []
)
),
),
(
"Annotation build",
_format_hint_values(
annotation_reference_summary.get(
"annotation_build_hints", []
)
),
),
(
"Reference provider",
_format_hint_values(
annotation_reference_summary.get(
"reference_provider_hints", []
)
),
),
(
"Annotation provider",
_format_hint_values(
annotation_reference_summary.get(
"annotation_provider_hints", []
)
),
),
]
DataTable.from_pandas(
pd.DataFrame(hint_rows, columns=["Evidence", "Hints"]),
paging=False,
searchable=False,
)
examples = annotation_summary.get("unstranded_examples") or []
if examples:
with h3("Unstranded Annotation Examples"):
pre("\n".join(examples))
# Add Bambu QC section with warnings
if bambu_qc:

370
bin/workflow_glue/tests/test_prepare_annotation_reference.py Normal file
View File

@ -0,0 +1,370 @@
"""Tests for annotation/reference preparation."""
import gzip
from pathlib import Path
from types import SimpleNamespace
import pytest
from workflow_glue import get_components
from workflow_glue import prepare_annotation_reference
def _write(path, text):
path.write_text(text, encoding="utf-8")
return path
def _write_gzip(path, text):
with gzip.open(path, "wt", encoding="utf-8") as handle:
handle.write(text)
return path
SIMPLE_REFERENCE = ">chr1\nAAAA\n"
SIMPLE_GTF = 'chr1\tsim\ttranscript\t1\t4\t.\t+\t.\tgene_id "g1"; transcript_id "t1";\n'
GFF3_INPUT = """##gff-version 3
chr1\tsim\tgene\t1\t10\t.\t+\t.\tID=gene1
chr1\tsim\tmRNA\t1\t10\t.\t+\t.\tID=transcript1;Parent=gene1
chr1\tsim\texon\t1\t5\t.\t+\t.\tParent=transcript1
chr1\tsim\texon\t6\t10\t.\t+\t.\tParent=transcript1
"""
GFF3_AS_GTF = (
'chr1\tsim\ttranscript\t1\t10\t.\t+\t.\tgene_id "gene1"; '
'transcript_id "transcript1";\n'
'chr1\tsim\texon\t1\t5\t.\t+\t.\tgene_id "gene1"; '
'transcript_id "transcript1";\n'
'chr1\tsim\texon\t6\t10\t.\t+\t.\tgene_id "gene1"; '
'transcript_id "transcript1";\n'
)
NCBI_REFERENCE = """>chr1 GRCh38 RefSeq primary assembly
AAAA
>chrExtra
TTTT
"""
NCBI_GTF = (
"#!genome-build GRCh38\n"
"#!genome-version GRCh38.p14 RefSeq\n"
'chr1\tRefSeq\ttranscript\t1\t4\t.\t-\t.\tgene_id "transcript_id"; '
'transcript_id "rna-XM_000001.1";\n'
'chr1\tRefSeq\texon\t1\t4\t.\t-\t.\ttranscript_id "rna-XM_000001.1"; '
'exon_number "1";\n'
)
MIXED_STRANDED_GTF = (
"# header\n"
'chr1\tsim\ttranscript\t1\t4\t.\t+\t.\tgene_id "g1"; transcript_id "t1";\n'
'chrMissing\tsim\ttranscript\t1\t4\t.\t.\t.\tgene_id "g2"; '
'transcript_id "t2";\n'
)
MOUSE_REFERENCE = ">chr1 dna:GRCm39 primary assembly\nAAAA\n"
MOUSE_GTF = (
"#!genome-build GRCm39\n"
'chr1\tGENCODE\ttranscript\t1\t4\t.\t+\t.\tgene_id "g1"; '
'transcript_id "t1";\n'
)
@pytest.mark.parametrize(
"ref_data,ref_name,annot_data,annot_name,gff_output",
[
(SIMPLE_REFERENCE, "reference.fa", SIMPLE_GTF,
"annotation.gtf", None),
(SIMPLE_REFERENCE, "reference.fa.gz", SIMPLE_GTF,
"annotation.gtf.gz", None),
(SIMPLE_REFERENCE, "reference.fa", GFF3_INPUT,
"annotation.gff3.gz", GFF3_AS_GTF),
(NCBI_REFERENCE, "reference.fna.gz", NCBI_GTF,
"annotation.gtf.gz", None),
(SIMPLE_REFERENCE + ">chrExtra\nTTTT\n", "reference.fa",
MIXED_STRANDED_GTF, "annotation.gtf", None),
(MOUSE_REFERENCE, "GRCm39.genome.fa", MOUSE_GTF,
"gencode.vM33.annotation.gtf", None),
],
ids=[
"simple_gtf", "gzipped_gtf", "gff3_input",
"ncbi_format", "mixed_stranded", "mouse_gencode"],
)
def test_prepare_all_formats_produce_valid_outputs(
ref_data, ref_name, annot_data, annot_name, gff_output, tmp_path, monkeypatch
):
"""All supported input formats produce annotation.gtf and reference.fasta."""
# write reference
ref_path = tmp_path / ref_name
if ref_name.endswith(".gz"):
_write_gzip(ref_path, ref_data)
else:
_write(ref_path, ref_data)
# write annotation
annot_path = tmp_path / annot_name
if annot_name.endswith(".gz"):
_write_gzip(annot_path, annot_data)
else:
_write(annot_path, annot_data)
# mock gffread (only invoked for GFF3 inputs)
if gff_output:
def fake_run(command, check, capture_output, text):
output_idx = command.index("-o") + 1
Path(command[output_idx]).write_text(gff_output, encoding="utf-8")
return SimpleNamespace(returncode=0, stdout="", stderr="")
monkeypatch.setattr(prepare_annotation_reference.subprocess, "run", fake_run)
out_dir = tmp_path / "prepared"
summary = prepare_annotation_reference.prepare_annotation_reference(
annot_path, ref_path, out_dir
)
# output files must exist
assert (out_dir / "annotation.gtf").exists()
assert (out_dir / "reference.fasta").exists()
assert (out_dir / "annotation_reference_summary.json").exists()
# paths returned in summary must match
assert summary["annotation"]["prepared"] == str(out_dir / "annotation.gtf")
assert summary["reference"]["prepared"] == str(out_dir / "reference.fasta")
# outputs must be valid for downstream tools (Bambu, SQANTI)
reference = Path(summary["reference"]["prepared"])
annotation = Path(summary["annotation"]["prepared"])
assert reference.read_text(encoding="utf-8").startswith(">")
assert summary["analysis_seqnames"]["has_overlap"] is True
# all annotation records must be stranded with gene_id and transcript_id
records = [
line.rstrip("\n").split("\t")
for line in annotation.read_text(encoding="utf-8").splitlines()
if line and not line.startswith("#")
]
assert records
for fields in records:
assert len(fields) == 9
assert fields[6] in {"+", "-"}
assert 'gene_id "' in fields[8]
assert 'transcript_id "' in fields[8]
def test_ncbi_format_sanitises_gene_ids(tmp_path):
"""NCBI gene_id='transcript_id' should be sanitised to actual transcript_id."""
reference = _write_gzip(tmp_path / "reference.fna.gz", NCBI_REFERENCE)
annotation = _write_gzip(tmp_path / "annotation.gtf.gz", NCBI_GTF)
summary = prepare_annotation_reference.prepare_annotation_reference(
annotation,
reference,
tmp_path / "prepared",
)
prepared_text = Path(summary["annotation"]["prepared"]).read_text(encoding="utf-8")
assert 'gene_id "transcript_id"' not in prepared_text
assert prepared_text.count('gene_id "rna-XM_000001.1"') == 2
assert summary["annotation"]["sanitised_attribute_records"] == 2
assert any(
"Sanitised 2 annotation records" in warning
for warning in summary["warnings"]
)
def test_gff3_conversion_via_gffread(tmp_path, monkeypatch):
"""GFF3 inputs invoke gffread and convert to GTF format."""
reference = _write(tmp_path / "reference.fa", SIMPLE_REFERENCE)
annotation = _write(tmp_path / "annotation.gff3", GFF3_INPUT)
captured = {}
def fake_run(command, check, capture_output, text):
captured["command"] = command
output_idx = command.index("-o") + 1
Path(command[output_idx]).write_text(GFF3_AS_GTF, encoding="utf-8")
return SimpleNamespace(returncode=0, stdout="", stderr="")
monkeypatch.setattr(prepare_annotation_reference.subprocess, "run", fake_run)
summary = prepare_annotation_reference.prepare_annotation_reference(
annotation,
reference,
tmp_path / "prepared",
gffread="custom-gffread",
)
assert summary["annotation"]["was_gff"] is True
assert captured["command"][0] == "custom-gffread"
assert "-T" in captured["command"]
prepared_text = Path(summary["annotation"]["prepared"]).read_text(encoding="utf-8")
assert 'gene_id "gene1"' in prepared_text
assert 'transcript_id "transcript1"' in prepared_text
assert prepared_text.count("exon") == 2
def test_unstranded_records_filtered_and_saved_separately(tmp_path):
"""Unstranded entries removed from main annotation, saved to separate file."""
reference = _write(tmp_path / "reference.fa", SIMPLE_REFERENCE)
annotation = _write(
tmp_path / "annotation.gtf",
SIMPLE_GTF
+ 'chr1\tsim\ttranscript\t5\t8\t.\t.\t.\t'
+ 'gene_id "g2"; transcript_id "t2";\n',
)
summary = prepare_annotation_reference.prepare_annotation_reference(
annotation,
reference,
tmp_path / "prepared",
)
prepared = Path(summary["annotation"]["prepared"]).read_text(encoding="utf-8")
unstranded = Path(
summary["annotation"]["unstranded_path"]).read_text(encoding="utf-8")
assert 'transcript_id "t1"' in prepared
assert 'transcript_id "t2"' not in prepared
assert 'transcript_id "t2"' in unstranded
assert summary["annotation"]["excluded_unstranded_records"] == 1
assert any("Excluded 1 unstranded" in w for w in summary["warnings"])
def test_seqname_warnings_reflect_filtered_annotation(tmp_path):
"""Seqname overlap checks use filtered annotation, not unfiltered."""
reference = _write(
tmp_path / "reference.fa", SIMPLE_REFERENCE + ">chrExtra\nTTTT\n")
annotation = _write(tmp_path / "annotation.gtf", MIXED_STRANDED_GTF)
summary = prepare_annotation_reference.prepare_annotation_reference(
annotation,
reference,
tmp_path / "prepared",
)
# chrMissing was filtered as unstranded, so NOT in seqname warnings
assert "chrMissing" not in summary["seqnames"]["only_in_annotation"]
assert "chrExtra" in summary["seqnames"]["only_in_reference"]
assert not any(
"chrMissing" in w and "seqnames" in w.lower()
for w in summary["warnings"]
)
def test_gffread_failure_raises_error(tmp_path, monkeypatch):
"""Failed GFF3 conversion should abort before downstream processing."""
reference = _write(tmp_path / "reference.fa", SIMPLE_REFERENCE)
annotation = _write(tmp_path / "annotation.gff3", GFF3_INPUT)
def fake_run(command, check, capture_output, text):
return SimpleNamespace(returncode=1, stdout="", stderr="gffread error")
monkeypatch.setattr(prepare_annotation_reference.subprocess, "run", fake_run)
with pytest.raises(ValueError, match="Failed to convert annotation"):
prepare_annotation_reference.prepare_annotation_reference(
annotation,
reference,
tmp_path / "prepared",
)
def test_no_seqname_overlap_raises_error(tmp_path):
"""Completely mismatched reference and annotation should abort."""
reference = _write(tmp_path / "reference.fa", ">chrOther\nAAAA\n")
annotation = _write(tmp_path / "annotation.gtf", SIMPLE_GTF)
with pytest.raises(ValueError, match="No overlapping seqnames"):
prepare_annotation_reference.prepare_annotation_reference(
annotation,
reference,
tmp_path / "prepared",
)
def test_malformed_gtf_raises_error(tmp_path):
"""GTF with wrong number of columns should abort."""
reference = _write(tmp_path / "reference.fa", SIMPLE_REFERENCE)
annotation = _write(
tmp_path / "annotation.gtf",
'chr1\tsim\ttranscript\t1\t4\t.\t+\t.\n',
)
with pytest.raises(ValueError, match="expected 9 GTF columns"):
prepare_annotation_reference.prepare_annotation_reference(
annotation,
reference,
tmp_path / "prepared",
)
def test_all_records_unstranded_raises_error(tmp_path):
"""If filtering removes all records, abort before downstream tools."""
reference = _write(tmp_path / "reference.fa", SIMPLE_REFERENCE)
annotation = _write(
tmp_path / "annotation.gtf",
'chr1\tsim\ttranscript\t1\t4\t.\t.\t.\tgene_id "g1"; transcript_id "t1";\n',
)
with pytest.raises(ValueError, match="reference annotation"):
prepare_annotation_reference.prepare_annotation_reference(
annotation,
reference,
tmp_path / "prepared",
)
def test_conflicting_build_hints_warn(tmp_path):
"""Mismatched genome build hints produce warning but do not fail."""
reference = _write(
tmp_path / "reference.fa",
">chr1 GRCh38 primary assembly\nAAAA\n",
)
annotation = _write(
tmp_path / "annotation.gtf",
"#!genome-build GRCh37\n" + SIMPLE_GTF,
)
summary = prepare_annotation_reference.prepare_annotation_reference(
annotation,
reference,
tmp_path / "prepared",
)
assert set(summary["reference_build_hints"]) == {"GRCh38"}
assert set(summary["annotation_build_hints"]) == {"GRCh37"}
assert any("different genome build hints" in w for w in summary["warnings"])
def test_conflicting_provider_hints_warn(tmp_path):
"""Mismatched provider hints produce warning but do not fail."""
reference = _write(tmp_path / "ensembl.GRCh38.fa", ">chr1\nAAAA\n")
annotation = _write(tmp_path / "refseq.GRCh38.annotation.gtf", SIMPLE_GTF)
summary = prepare_annotation_reference.prepare_annotation_reference(
annotation,
reference,
tmp_path / "prepared",
)
assert any("different provider hints" in w for w in summary["warnings"])
def test_filename_hints_captured(tmp_path):
"""Build and provider hints extracted from filenames."""
reference = _write(tmp_path / "reference.hg38.fa", ">chr1\nAAAA\n")
annotation = _write(tmp_path / "gencode.GRCh37.annotation.gtf", SIMPLE_GTF)
summary = prepare_annotation_reference.prepare_annotation_reference(
annotation,
reference,
tmp_path / "prepared",
)
assert "GRCh38" in summary["reference_build_hints"]
assert "GRCh37" in summary["annotation_build_hints"]
assert "GENCODE" in summary["annotation_provider_hints"]
def test_cli_component_discoverable():
"""workflow-glue discovers prepare_annotation_reference command."""
components = get_components(allowed_components=["prepare_annotation_reference"])
assert "prepare_annotation_reference" in components

30
bin/workflow_glue/tests/test_report.py
View File

@ -44,15 +44,17 @@ def _write(path, text):
def test_report_main_accepts_optional_file_sentinels(monkeypatch, tmp_path):
"""The report entry point should tolerate null-object sentinel files."""
tables = []
monkeypatch.setattr(report.labs, "LabsReport", _FakeReport)
monkeypatch.setattr(report, "Tabs", _FakeTabs)
monkeypatch.setattr(report, "p", lambda *args, **kwargs: None)
monkeypatch.setattr(report, "pre", lambda *args, **kwargs: None)
monkeypatch.setattr(report, "_create_warning_banner", lambda *args, **kwargs: None)
monkeypatch.setattr(report.fastcat, "SeqSummary", lambda *args, **kwargs: None)
monkeypatch.setattr(
report.DataTable,
"from_pandas",
staticmethod(lambda *args, **kwargs: None),
staticmethod(lambda table, *args, **kwargs: tables.append(table.copy())),
)
metadata = _write(
@ -66,6 +68,27 @@ def test_report_main_accepts_optional_file_sentinels(monkeypatch, tmp_path):
cohort = tmp_path / "cohort"
cohort.mkdir()
reference = cohort / "reference"
reference.mkdir()
_write(
reference / "annotation_reference_summary.json",
json.dumps({
"seqname_overlap": ["chr1"],
"only_in_annotation": ["chrMissing"],
"only_in_reference": ["chrExtra"],
"annotation": {
"kept_records": 10,
"excluded_unstranded_records": 2,
"sanitised_attribute_records": 1,
"unstranded_examples": ["chr1\tsim\ttranscript\t1\t4\t.\t.\t."],
},
"reference_build_hints": ["GRCh38"],
"annotation_build_hints": [],
"reference_provider_hints": [],
"annotation_provider_hints": [],
"warnings": ["Warning: Some seqnames are present in the annotation."],
}),
)
samples = tmp_path / "samples"
samples.mkdir()
@ -103,6 +126,11 @@ def test_report_main_accepts_optional_file_sentinels(monkeypatch, tmp_path):
report.main(args)
assert out_report.exists()
assert any("Overlapping seqnames" in table.to_string() for table in tables)
assert any(
"Annotation attributes sanitised" in table.to_string() for table in tables
)
assert any("GRCh38" in table.to_string() for table in tables)
def test_pychopper_tables_uses_sample_directory_names(tmp_path):

42
bin/workflow_glue/tests/test_report_qc.py
View File

@ -77,6 +77,48 @@ def test_load_de_qc_missing(tmp_path):
assert result is None
def test_load_annotation_reference_summary_exists(tmp_path):
"""Test loading reference/annotation prep summary when file exists."""
from workflow_glue.report import _load_annotation_reference_summary
summary_dir = tmp_path / "cohort" / "reference"
summary_dir.mkdir(parents=True)
summary_data = {
"seqname_overlap": ["chr1"],
"only_in_annotation": ["chrMissing"],
"only_in_reference": ["chrExtra"],
"warnings": ["Warning: test"],
}
with open(summary_dir / "annotation_reference_summary.json", "w") as f:
json.dump(summary_data, f)
result = _load_annotation_reference_summary(tmp_path / "cohort")
assert result is not None
assert result["seqname_overlap"] == ["chr1"]
assert result["only_in_annotation"] == ["chrMissing"]
def test_load_annotation_reference_summary_missing(tmp_path):
"""Test loading reference/annotation prep summary when file is missing."""
from workflow_glue.report import _load_annotation_reference_summary
cohort_dir = tmp_path / "cohort"
cohort_dir.mkdir()
result = _load_annotation_reference_summary(cohort_dir)
assert result is None
def test_format_hint_values():
"""Build/provider hints should be compactly formatted for the report."""
from workflow_glue.report import _format_hint_values
assert _format_hint_values([]) == "None detected"
assert _format_hint_values(["GRCh38"]) == "GRCh38"
assert _format_hint_values(["GRCh38", "GENCODE"]) == "GRCh38, GENCODE"
def test_warning_banner_creation():
"""Test that warning banner can be created."""
from workflow_glue.report import _create_warning_banner

2
docs/08_outputs.md
View File

@ -10,6 +10,8 @@ Output files may be aggregated including information for all samples or provided
| Aligned BAM | cohort/alignments/{{ alias }}.aligned.sorted.bam | Genome-aligned BAM used for bambu, optional SQANTI3 QC, and IGV. | per-sample |
| Aligned BAM index | cohort/alignments/{{ alias }}.aligned.sorted.bam.bai | Index for the aligned BAM. | per-sample |
| Alignment summary | cohort/alignments/{{ alias }}.flagstat.txt | samtools flagstat output for the aligned BAM. | per-sample |
| Reference and annotation preparation summary | cohort/reference/annotation_reference_summary.json | Summary of reference and annotation preparation, including seqname overlap, build/provider hints, and excluded unstranded annotation counts. | aggregated |
| Excluded unstranded annotation records | cohort/reference/unstranded_annotation.gtf | Full set of annotation records excluded because their strand was not '+' or '-'. Present only when unstranded records are found. | aggregated |
| Cohort transcriptome GTF | cohort/transcripts.gtf | Joint bambu transcript model used as the primary cohort transcriptome. | aggregated |
| Cohort transcriptome FASTA | cohort/cohort.transcriptome.fa | Transcript sequences derived from the joint cohort GTF. | aggregated |
| Cohort transcript counts | cohort/transcript_counts.tsv | Transcript-level count matrix produced by bambu. | aggregated |

2
main.nf
View File

@ -275,6 +275,8 @@ workflow pipeline {
.concat(report.report.map { [it, null] })
.concat(workflow_params.map { [it, null] })
.concat(pychopper_dirs.map { meta, pychopper_dir -> [pychopper_dir, "ingress_results/${meta.alias}"] })
.concat(transcriptome.annotation_reference_summary.map { [it, "cohort/reference"] })
.concat(transcriptome.unstranded_annotation.map { [it, "cohort/reference"] })
.concat(transcriptome.joint_gtf.map { [it, "cohort"] })
.concat(transcriptome.joint_fasta.map { [it, "cohort"] })
.concat(transcriptome.joint_transcript_counts.map { [it, "cohort"] })

16
output_definition.json
View File

@ -64,6 +64,22 @@
"optional": false,
"type": "per-sample"
},
"annotation-reference-summary": {
"filepath": "cohort/reference/annotation_reference_summary.json",
"title": "Reference and annotation preparation summary",
"description": "Summary of reference and annotation preparation, including seqname overlap, build/provider hints, and excluded unstranded annotation counts.",
"mime-type": "application/json",
"optional": false,
"type": "aggregated"
},
"unstranded-annotation": {
"filepath": "cohort/reference/unstranded_annotation.gtf",
"title": "Excluded unstranded annotation records",
"description": "Full set of annotation records excluded because their strand was not '+' or '-'. Present only when unstranded records are found.",
"mime-type": "text/plain",
"optional": true,
"type": "aggregated"
},
"cohort-gtf": {
"filepath": "cohort/transcripts.gtf",
"title": "Cohort transcriptome GTF",

158
subworkflows/transcriptome.nf
View File

@ -3,114 +3,26 @@ nextflow.enable.dsl = 2
OPTIONAL_FILE = file("$projectDir/data/OPTIONAL_FILE")
process decompressReference {
label "wf_transcriptomes"
cpus 1
memory "2 GB"
input:
path compressed_ref
output:
path "${compressed_ref.baseName}", emit: decompressed_ref
script:
"""
gzip -dc "${compressed_ref}" > "${compressed_ref.baseName}"
"""
}
process decompressAnnotation {
label "wf_transcriptomes"
cpus 1
memory "2 GB"
input:
path compressed_annotation
output:
path "${compressed_annotation.baseName}", emit: decompressed_annotation
script:
"""
gzip -dc "${compressed_annotation}" > "${compressed_annotation.baseName}"
"""
}
process normaliseAnnotationToGtf {
process prepareAnnotationReference {
label "wf_transcriptomes"
cpus 1
memory "4 GB"
input:
path annotation
output:
path "annotation.gtf", emit: gtf
script:
"""
gffread -T "${annotation}" -o annotation.gtf
"""
}
process filterUnstrandedAnnotation {
label "wf_transcriptomes"
cpus 1
memory "2 GB"
input:
path "annotation.gtf"
output:
tuple stdout, path("annotation_stranded.gtf"), emit: filtered
script:
"""
awk '
BEGIN { OFS = "\\t" }
/^#/ { print; next }
(\$7 == "+" || \$7 == "-") { print; next }
{ print \$0 >> "unstranded.gtf" }
' annotation.gtf > annotation_stranded.gtf
if [ -s unstranded.gtf ]; then
echo "Warning: Unstranded entries found and excluded from differential expression analysis."
echo "If running with reference-guided transcriptome source, consider increasing read depth to reduce unstranded annotations."
echo "If running with precomputed transcriptome source, ensure your ref_annotation gtf file contains only '+' or '-' strand entries."
echo "A sample of unstranded entries:"
head -n 20 unstranded.gtf
fi
"""
}
process validateReferenceAnnotation {
label "wf_transcriptomes"
cpus 1
memory "2 GB"
input:
path "annotation.gtf"
path "reference.fasta"
path ref_annotation
path ref_genome
output:
stdout emit: warnings
path "annotation.gtf", emit: annotation
path "reference.fasta", emit: reference
path "annotation_reference_summary.json", emit: summary
path "unstranded_annotation.gtf", optional: true, emit: unstranded
script:
"""
grep -v '^#' annotation.gtf | cut -f1 | sort -u > annotation_ids.txt
awk '/^>/ {print substr(\$1,2)}' reference.fasta | sort -u > reference_ids.txt
matches=\$(comm -12 annotation_ids.txt reference_ids.txt || true)
only_in_annotation=\$(comm -23 annotation_ids.txt reference_ids.txt || true)
only_in_reference=\$(comm -13 annotation_ids.txt reference_ids.txt || true)
if [[ -z "\$matches" ]]; then
echo "ERROR: No overlapping seqnames were found between the reference annotation and the reference genome." >&2
echo "Annotation ID examples:" >&2
head -n 5 annotation_ids.txt >&2
echo "Reference ID examples:" >&2
head -n 5 reference_ids.txt >&2
exit 78
fi
if [[ -n "\$only_in_annotation" ]]; then
echo "Warning: Some seqnames are present in the annotation but not the genome:"
echo "\$only_in_annotation" | head -n 5
fi
if [[ -n "\$only_in_reference" ]]; then
echo "Warning: Some seqnames are present in the genome but not the annotation:"
echo "\$only_in_reference" | head -n 5
fi
workflow-glue prepare_annotation_reference \
--annotation "${ref_annotation}" \
--reference "${ref_genome}" \
--out_dir prepared
mv prepared/* .
"""
}
@ -272,14 +184,9 @@ process runJointSqanti {
String skip_orf = params.sqanti_skip_orf ? "--skipORF" : ""
"""
mkdir sqanti_cohort
awk '
BEGIN { OFS = "\\t" }
/^#/ { print; next }
(\$7 == "+" || \$7 == "-") { print; next }
' "${annotation}" > reference_annotation_stranded.gtf
sqanti3_qc.py \
--isoforms "${gtf}" \
--refGTF "reference_annotation_stranded.gtf" \
--refGTF "${annotation}" \
--refFasta "${reference}" \
${skip_orf} \
--force_id_ignore \
@ -310,14 +217,9 @@ process runPerSampleSqanti {
String skip_orf = params.sqanti_skip_orf ? "--skipORF" : ""
"""
mkdir "${meta.alias}_sqanti"
awk '
BEGIN { OFS = "\\t" }
/^#/ { print; next }
(\$7 == "+" || \$7 == "-") { print; next }
' "${annotation}" > reference_annotation_stranded.gtf
sqanti3_qc.py \
--isoforms "${gtf}" \
--refGTF "reference_annotation_stranded.gtf" \
--refGTF "${annotation}" \
--refFasta "${reference}" \
${skip_orf} \
--force_id_ignore \
@ -370,34 +272,14 @@ workflow transcriptome_analysis {
ref_annotation
sample_sheet
main:
analysis_reference = ref_genome
if (params.ref_genome.toLowerCase().endsWith("gz")) {
decompressed_reference = decompressReference(ref_genome)
analysis_reference = decompressed_reference.decompressed_ref
}
analysis_annotation = ref_annotation
if (params.ref_annotation.toLowerCase().endsWith("gz")) {
decompressed_annotation = decompressAnnotation(ref_annotation)
analysis_annotation = decompressed_annotation.decompressed_annotation
}
if (params.ref_annotation.toLowerCase() ==~ /.*\.gff3?(\.gz)?$/) {
normalised_annotation = normaliseAnnotationToGtf(analysis_annotation)
analysis_annotation = normalised_annotation.gtf
}
filtered_annotation = filterUnstrandedAnnotation(analysis_annotation)
analysis_annotation = filtered_annotation.filtered.map { warning_text, annotation ->
if (warning_text?.trim()) {
log.warn(warning_text.trim())
}
annotation
}
validateReferenceAnnotation(analysis_annotation, analysis_reference).map { stdoutput ->
prepared_reference_annotation = prepareAnnotationReference(ref_annotation, ref_genome)
prepared_reference_annotation.warnings.map { stdoutput ->
if (stdoutput) {
log.warn(stdoutput)
log.warn(stdoutput.trim())
}
}
analysis_annotation = prepared_reference_annotation.annotation
analysis_reference = prepared_reference_annotation.reference
genome_index = buildMinimapIndex(analysis_reference)
@ -442,6 +324,8 @@ workflow transcriptome_analysis {
emit:
reference = ref_genome
annotation = ref_annotation
annotation_reference_summary = prepared_reference_annotation.summary
unstranded_annotation = prepared_reference_annotation.unstranded
alignments = aligned.bam
joint_dir = joint_bambu.dir
joint_gtf = joint_bambu.gtf