diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index eb48389..d3e38c9 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -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" diff --git a/README.md b/README.md index 5af3275..67fcf3e 100644 --- a/README.md +++ b/README.md @@ -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 | diff --git a/bin/workflow_glue/prepare_annotation_reference.py b/bin/workflow_glue/prepare_annotation_reference.py new file mode 100644 index 0000000..efab524 --- /dev/null +++ b/bin/workflow_glue/prepare_annotation_reference.py @@ -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"(?"): + 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 diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py index e77a126..d089f6e 100644 --- a/bin/workflow_glue/report.py +++ b/bin/workflow_glue/report.py @@ -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: diff --git a/bin/workflow_glue/tests/test_prepare_annotation_reference.py b/bin/workflow_glue/tests/test_prepare_annotation_reference.py new file mode 100644 index 0000000..7b29f69 --- /dev/null +++ b/bin/workflow_glue/tests/test_prepare_annotation_reference.py @@ -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 diff --git a/bin/workflow_glue/tests/test_report.py b/bin/workflow_glue/tests/test_report.py index afddd9b..8c58242 100644 --- a/bin/workflow_glue/tests/test_report.py +++ b/bin/workflow_glue/tests/test_report.py @@ -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): diff --git a/bin/workflow_glue/tests/test_report_qc.py b/bin/workflow_glue/tests/test_report_qc.py index 134d2ce..a9dd6cf 100644 --- a/bin/workflow_glue/tests/test_report_qc.py +++ b/bin/workflow_glue/tests/test_report_qc.py @@ -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 diff --git a/docs/08_outputs.md b/docs/08_outputs.md index 041f683..bf9de17 100644 --- a/docs/08_outputs.md +++ b/docs/08_outputs.md @@ -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 | diff --git a/main.nf b/main.nf index c3e679c..5e12a10 100644 --- a/main.nf +++ b/main.nf @@ -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"] }) diff --git a/output_definition.json b/output_definition.json index 4cf859f..88895e4 100644 --- a/output_definition.json +++ b/output_definition.json @@ -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", diff --git a/subworkflows/transcriptome.nf b/subworkflows/transcriptome.nf index 2ca14bd..f473b4d 100644 --- a/subworkflows/transcriptome.nf +++ b/subworkflows/transcriptome.nf @@ -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