diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index a81a490..8550b90 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -4,10 +4,9 @@ include: file: "wf-containers.yaml" variables: - NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -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 -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_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/hg38_chr20.fa \ - --transcriptome_source reference-guided \ --ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.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" @@ -22,13 +21,14 @@ macos-run: - x86 -# run all our test in AWS, cos EKS runners don't work docker-run: - rules: - - when: never - -aws-run: - stage: test + tags: + - linux + - prod + - amd64 + - eks + - xlarge-highio + - docker artifacts: when: always paths: @@ -59,7 +59,10 @@ aws-run: "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" + "mismatch_seqid_error","mismatch_seqid_warn", + "codex_discover", "codex_fixed_annotation", "codex_direct_rna", + "codex_de_analysis", "codex_invalid_mode", + "codex_conflicting_library_flags" ] rules: # NOTE As we're overriding the rules block for the included docker-run @@ -69,44 +72,45 @@ aws-run: when: never - if: $MATRIX_NAME == "isoforms" variables: - NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -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 -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 --transcriptome_source reference-guided \ - --ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa --ref_annotation ${CI_PROJECT_NAME}/data/chr20/gencode.v22.annotation.chr20.gtf --pychopper_backend phmm" + 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 -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 -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 --transcriptome_source reference-guided \ - --ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa --ref_annotation ${CI_PROJECT_NAME}/data/chr20/gencode.v22.annotation.chr20.gtf --pychopper_backend phmm \ + 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" variables: - NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -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 -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 --transcriptome_source reference-guided \ + 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 \ -c ${CI_PROJECT_NAME}/data/demo.nextflow.config " - NF_IGNORE_PROCESSES: run_gffcompare,filter_unstranded_annotation,validate_ref_annotation,preprocess_reads,faidx,gz_faidx,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome + NF_IGNORE_PROCESSES: preprocess_reads + ASSERT_NEXTFLOW_FAILURE: "1" + ASSERT_NEXTFLOW_FAILURE_REXP: "Missing required parameter: --ref_annotation" - if: $MATRIX_NAME == "differential_expression" variables: - NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -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 -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_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/hg38_chr20.fa --transcriptome_source reference-guided \ + --ref_genome ${CI_PROJECT_NAME}/data/differential_expression/hg38_chr20.fa \ --ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.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 - if: $MATRIX_NAME == "only_differential_expression" variables: - NF_BEFORE_SCRIPT: "mkdir -p ${CI_PROJECT_NAME}/data/ && wget -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 -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_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_source precomputed \ + --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' \ - --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression/ref_transcriptome.fasta \ --sample_sheet test_data/sample_sheet.csv \ -c ${CI_PROJECT_NAME}/data/demo.nextflow.config " NF_IGNORE_PROCESSES: > @@ -114,14 +118,13 @@ aws-run: 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 -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 -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_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_source precomputed \ + --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' \ - --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression/ref_transcriptome.fasta \ --sample_sheet test_data/sample_sheet.csv \ -c ${CI_PROJECT_NAME}/data/demo.nextflow.config " NF_IGNORE_PROCESSES: > @@ -129,7 +132,7 @@ aws-run: 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 -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 -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_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 \ @@ -140,30 +143,30 @@ aws-run: -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 + 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 -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 -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_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_source precomputed --de_analysis \ + --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 --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_rna.fna.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 + 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 -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 -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_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 --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression/Homo_sapiens.GRCh38.cdna.all.fa.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: > @@ -171,24 +174,24 @@ aws-run: 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 -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 -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_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_source precomputed --de_analysis \ + --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 --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression_mouse/gencode.vM33.transcripts.fa.gz \ + --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 + 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 -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 -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_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 --transcriptome_source reference-guided \ + --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 " @@ -196,13 +199,12 @@ aws-run: 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 -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 -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_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/hg38_chr20.fa \ --ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.gff \ - --direct_rna --minimap2_index_opts '-k 15' \ - --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression/ref_transcriptome.fasta \ + --direct_rna --transcriptome_mode fixed_annotation --minimap2_index_opts '-k 15' \ --sample_sheet test_data/sample_sheet.csv \ --igv \ -c ${CI_PROJECT_NAME}/data/demo.nextflow.config " @@ -211,38 +213,38 @@ aws-run: 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 -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 -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_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_source precomputed --de_analysis \ + --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 --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_rna.fna.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 -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 -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: Seqid mismatch found between the provided ref_annotation + 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 -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 -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_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 \ @@ -253,7 +255,51 @@ aws-run: -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 'are only present in the reference genome' '.nextflow.log'" - - + 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" + 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" + NF_IGNORE_PROCESSES: preprocess_reads + AFTER_NEXTFLOW_CMD: > + 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" + 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" + NF_IGNORE_PROCESSES: preprocess_reads + 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" + 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" + NF_IGNORE_PROCESSES: preprocess_reads + 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" + 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" + NF_IGNORE_PROCESSES: preprocess_reads + 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" + 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" + 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" + ASSERT_NEXTFLOW_FAILURE: "1" + ASSERT_NEXTFLOW_FAILURE_REXP: ".*--cdna_preprocess cannot be used together with --direct_rna" diff --git a/CHANGELOG.md b/CHANGELOG.md index 29f9d61..d160211 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -99,10 +99,9 @@ This patch release of wf-transcriptomes updates internal workflow naming, and do ## [v1.2.0] ### Added - Workflow now accepts BAM or FASTQ files as input (using the --bam or --fastq parameters, respectively). +- Error message when running in `de_analysis` mode and `ref_annotation` input file contains unstranded annotations. ### Changed - MA plot in the `results_dge.pdf` has been updated to match the MA plot in the report. -### Added -- Error message when running in `de_analysis` mode and `ref_annotation` input file contains unstranded annotations. ## [v1.1.1] ### Changed @@ -119,6 +118,8 @@ This patch release of wf-transcriptomes updates internal workflow naming, and do - `de_analysis/all_counts_filtered.tsv` to `de_analysis/filtered_transcript_counts_with_genes.tsv` - `de_analysis/de_tpm_transcript_counts.tsv` to `de_analysis/unfiltered_tpm_transcript_counts.tsv` - Minimum memory requirements to `32 GB`. +- More memory assigned to the fusion detection process. +- When no `--ref_annotation` is provided the workflow will still run but the output transcripts will not be annotated. However `--de_analysis` mode still requires a `--ref_annotation`. ### Added - Published isoforms table to output directory. - Output additional `de_analysis/cpm_gene_counts.tsv` with counts per million gene counts. @@ -126,9 +127,6 @@ This patch release of wf-transcriptomes updates internal workflow naming, and do - Add gene name column to the de_analysis counts TSV files. ### Fixed - Mapping stage using a single thread only. -### Changed -- More memory assigned to the fusion detection process. -- When no `--ref_annotation` is provided the workflow will still run but the output transcripts will not be annotated. However `--de_analysis` mode still requires a `--ref_annotation`. ## [v1.0.0] ### Added diff --git a/README.md b/README.md index 0d0664e..5af3275 100644 --- a/README.md +++ b/README.md @@ -1,18 +1,42 @@ # Transcriptomes -Transcriptome analysis of cDNA and direct RNA sequencing data. +Long-read transcriptome analysis using bambu with optional SQANTI3 QC, DESeq2, and DEXSeq. ## Introduction -This workflow can be used for the following: +This workflow analyses Oxford Nanopore long-read RNA sequencing data. It uses +[`bambu`](https://bioconductor.org/packages/bambu/) to build and quantify +transcript models, can optionally run +[`SQANTI3`](https://github.com/ConesaLab/SQANTI3) for transcript classification +and QC, and can optionally run +[`DESeq2`](https://bioconductor.org/packages/DESeq2/) and +[`DEXSeq`](https://bioconductor.org/packages/DEXSeq/) for differential +analysis. + +The workflow supports: + ++ transcript identification from either cDNA or direct RNA reads ++ optional cDNA preprocessing with `pychopper` before alignment ++ transcript discovery guided by a supplied genome and annotation ++ quantification against a supplied reference annotation ++ optional transcript classification and QC with `SQANTI3` ++ differential gene expression with `DESeq2` ++ differential transcript usage with `DEXSeq` + +The main transcriptome result is a shared `bambu` model built from all samples +together. The workflow also produces separate per-sample transcriptomes, so +each sample has its own GTF, FASTA, count tables, and optional `SQANTI3` +summary alongside the shared results. + +For users familiar with earlier transcriptome workflows, the main change is +that transcript discovery, quantification, and optional differential analysis +now use the shared `bambu` outputs rather than the older +StringTie/GffCompare/Salmon-based approach. The rest of this README explains +the current workflow in plain terms, while the `FAQ` and `Troubleshooting` +sections call out the main differences from the previous workflow version. -+ Identify RNA transcripts using either cDNA or direct RNA reads. -+ Reference aided transcriptome assembly. -+ Annotation of assembled transcripts. -+ Differential gene expression analysis using a pre-computed or assembled reference transcriptome. -+ Differential transcript usage analysis using a precomputed or assembled reference transcriptome. @@ -28,7 +52,7 @@ Minimum requirements: + CPUs = 8 + Memory = 32GB -Approximate run time: 15 minutes per sample, with 1 million reads and recommended resources. +Approximate run time: Varies with read depth and sample count; expect a small single-sample run to finish in under 30 minutes with the recommended resources. ARM processor support: False @@ -117,85 +141,155 @@ input_reads.fastq ─── input_directory ─── input_directory │ ├── reads1.fastq │ └── reads2.fastq └── barcode03 - └── reads0.fastq + └── reads0.fastq ``` + ## Pipeline overview -### 1. Concatenate input files and generate per read stats. -The [fastcat](https://github.com/epi2me-labs/fastcat) tool is used to concatenate multifile samples to be processed by the workflow. It will also output per read stats including average read lengths and qualities. +### 0. Background. -### 2. Preprocess cDNA. -If input sequences are cDNA [Pychopper](https://github.com/epi2me-labs/pychopper) is used to orient, trim and rescue full length cDNA reads and associated statistics. If the `direct_rna` parameter is selected this step will be skipped. +The methodology implemented within the wf-transcriptomics workflow follows from the largest independent long-read RNA benchmark to date. +The [Systematic assessment of long-read RNA-seq methods for transcript identification and quantification](https://www.nature.com/articles/s41592-024-02298-3) concluded that, in well-annotated genomes, reference-based methods perform best. +Our own previous research, benchmarking, and support of community members has shown that an automated, hands-off de-novo discovery pipeline to be bothersome for many use cases. +The wf-transcriptomes workflow therefore focuses on a reference-guided approach rather than a novelty-first one. -### 3. Build transcriptome. -If the `transcriptome_source` parameter is "reference-guided" a transcriptome will be built for each sample as outlined below. If the `transcriptome_source` is "precomputed" and the `reference_transcriptome` parameter is provided the workflow will skip step 3. +The benchmark paper above explicitly recommends `bambu` for identifying sample-specific transcriptomes in well-annotated organisms when only limited novelty is expected. +The paper also names `bambu` as one of the best options when quantification is important, which supports using it as the core engine for downstream DGE and DTU analyses. -#### 3.1 Align reads with reference genome. -The reference genome will be indexed and aligned using [Minimap2](https://github.com/lh3/minimap2). The output is sorted and converted to a BAM file using [Samtools](https://www.htslib.org/). Alignment stats are created from these using [Seqkit BAM](https://bioinf.shenwei.me/seqkit/usage/#bam). +In spike-in evaluations, `bambu` generally showed high precision and was among the better F1 performers. +This is an acceptable tradeoff for a production workflow where false transcript calls might be confound downstream analysis. +Users interested more in novel discovery may wish to amend the parameters of the workflow away from their defaults. +`bambu` also performed especially well on long non-spliced SIRVs, which supports its use on long-read datasets where transcript-end definition matters. -Additionally, the workflow will generate an IGV configuration file if `--igv` is selected. This file allows the user to view the aligned BAM in the EPI2ME Desktop Application in the Viewer tab. +The workflow's choice of SQANTI3 as a companion QC and annotation layer matches the benchmark paper, which used SQANTI3 categories and metrics as its transcript assessment framework; so our outputs align with the field’s standard reporting language. -#### 3.2 Chunk BAM -The aligned BAMs are split into chunks using the bundle_min_reads parameter (default: 50000). -#### 3.3 Assemble transcripts -[StringTie](https://ccb.jhu.edu/software/stringtie/) is then used to assemble the transcripts using the aligned segments in the chunked BAM files. The assembled transcript will be output as a [GFF file](https://www.ensembl.org/info/website/upload/gff3.html). If a `ref_annotation` file is provided this will also be included in the GFF. -#### 3.4 Merge Chunks -Transcript GFF files from the chunks with the same sample aliases will then be merged. +### 1. Getting your files into the workflow -#### 3.5 Annnotate -[GffCompare](https://ccb.jhu.edu/software/stringtie/gffcompare.html) is then used to compare query and reference annotations, merging records where appropriate and then annotating them. This also creates estimates of accuracy of the GFF files output in a stats file per sample. +The shared EPI2ME input handling collects FASTQ or BAM inputs, works out +whether you have a single sample or a multiplexed run, and produces per-sample +FASTQ files plus read statistics. These files are published under +`ingress_results//` and are used in the downstream report. -#### 3.6 Create transcriptomes -[Gffread](https://github.com/gpertea/gffread) is used to create a transcriptome FASTA file from the final GFF as well as a merged transcriptome that includes annotations in the FASTA headers where available. +### 2. Sample sheet formulation -### 4. Differential expression analysis +The sample sheet is optional for simple single-sample runs, but it becomes the +main source of sample names for multiplexed runs and is required for +`--de_analysis`. -Differential gene expression (DGE) and differential transcript usage (DTU) analyses aim to identify genes and transcripts that show statistically altered expression patterns. ++ Every row must contain `barcode` and `alias`. ++ `barcode` must use the usual ONT-style naming such as `barcode01`, + `barcode02`, and the values must be unique. ++ `alias` is the user-facing sample name, must be unique, and must not begin + with the word `barcode`. ++ If a `type` column is present, it must use one of: + `test_sample`, `positive_control`, `negative_control`, or + `no_template_control`. ++ If an `analysis_group` column is present, every row must have a value. ++ For `--de_analysis`, the sheet must also contain the primary condition + column, `condition` by default, plus any columns named in `--covariates`. -Differential Expression requires at least 2 replicates of each sample to compare (but we recommend three). You can see an example sample_sheet.csv below. +When multiplexed input folders are named by barcode, the workflow matches those +folder names against the `barcode` column. If the folders are named by alias, +the workflow can match them against `alias`, but the sample sheet still needs a +`barcode` column because the shared validator expects it. -#### Sample sheet condition column -The sample sheet should be a comma separated values file (.csv) and include at least three columns named `barcode`, `alias` and `condition`. -- Each `barcode` should refer to a directory of the same name in the input FASTQ directory (in the example below `barcode01` to `barcode06` reflect the `test_data` directory). -- The `alias` column allows you to rename each barcode to an alias that will be used in the report and other output files. -- The condition column will need to contain one of two keys to indicate the two samples being compared. Control must be one of the keys, used to indicate which samples will be used as the reference in the differential expression analysis. +Example sample sheet: -eg. sample_sheet.csv -``` -barcode,alias,condition -barcode01,sample01,control -barcode02,sample02,control -barcode03,sample03,control -barcode04,sample04,treated -barcode05,sample05,treated -barcode06,sample06,treated +```csv +barcode,alias,type,condition,batch +barcode01,control_rep1,test_sample,control,b1 +barcode02,control_rep2,test_sample,control,b2 +barcode03,treated_rep1,test_sample,treated,b1 +barcode04,treated_rep2,test_sample,treated,b2 ``` -#### 4.1 Merge cross sample transcriptomes -If a `ref_transcriptome` is not provided, the transcriptomes created by the workflow will be used for DE analysis. To do this, the GFF outputs of GffCompare are merged using StringTie. A final non redundant FASTA file of the transcripts is created using the merged GFF file and the reference genome using seqkit. +This example is suitable for a multiplexed run and also satisfies the minimum +requirements for a two-group DE/DTU comparison. -#### 4.2 Create a final non redundant transcriptome -The reads from all the samples will be aligned with the final non redundant transcriptome using Minimap2 in a splice aware manner. +### 3. Optional cDNA preprocessing -#### 4.3 Count genes and transcripts -[Salmon](https://github.com/COMBINE-lab/salmon) is used for transcript quantification, giving gene and transcript counts. +When `--cdna_preprocess` is enabled for cDNA libraries, the workflow runs +`pychopper` before alignment to classify, orient, and trim full-length reads. +This preprocessing stage is controlled by `--cdna_kit`, +`--pychopper_backend`, and optional extra `--pychopper_opts`, and its outputs +are published alongside the ingress results for each sample. -#### 4.4 edgeR based differential expression analysis -A statistical analysis is first performed using [edgeR](https://bioconductor.org/packages/release/bioc/html/edgeR.html) to identify the subset of differentially expressed genes using the gene counts as input. A normalisation factor is calculated for each sequence library using the default TMM method (see [McCarthy et al. (2012)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3378882/) for further details). The defined experimental design is used to calculate estimates of dispersion for each of the gene features. Statistical tests are calculated using the contrasts defined in the experimental design. The differentially expressed genes are corrected for false discovery (FDR) using the method of Benjamini & Hochberg ([Benjamini and Hochberg (1995)](https://www.jstor.org/stable/2346101)) +### 4. Genome alignment -#### 4.5 Pre-filtering of quantitative data using DRIMSeq -[DRIMSeq](https://bioconductor.org/packages/release/bioc/html/DRIMSeq.html) is used to filter the transcript count data from the Salmon analysis for differential transcript usage (DTU) analysis. The filter step will be used to select for genes and transcripts that satisfy rules for the number of samples in which a gene or transcript must be observed, and minimum threshold levels for the number of observed reads. The parameters used for filtering are `min_samps_gene_expr`, `min_samps_feature_expr`, `min_gene_expr`, and `min_feature_expr`. By default, any transcripts with zero expression or one transcript in all samples are filtered out at this stage. +Each sample is aligned to the supplied reference genome with +[`minimap2`](https://github.com/lh3/minimap2), then sorted and indexed with +[`samtools`](https://www.htslib.org/). The aligned BAMs under +`cohort/alignments/` are the main alignment files used for transcriptome +analysis, optional `SQANTI3` QC, and optional IGV viewing. -#### 4.6 Differential transcript usage using DEXSeq -Differential transcript usage analysis is performed using the R [DEXSeq](https://bioconductor.org/packages/release/bioc/html/DEXSeq.html) package ([Anders et al. (2012)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3460195/)). Similar to the edgeR package, DEXSeq estimates the variance between the biological replicates and applies generalised linear models for the statistical testing. The key difference is that the DEXSeq method looks for differences at the exon count level. DEXSeq uses the filtered transcript count data prepared earlier in this analysis. +### 5. Cohort transcriptome construction -#### 4.7 StageR stage-wise analysis of DGE and DTU -The final component of this isoform analysis is a stage-wise statistical test using the R software package [stageR](https://bioconductor.org/packages/release/bioc/html/stageR.html)([Van den Berge and Clement (2018)](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1277-0)). stageR uses (1) the raw p-values for DTU from the DEXSeq analysis in the previous section and (2) a false-discovery corrected set of p-values from testing whether individual genes contain at least one exon showing DTU. A hierarchical two-stage statistical testing evaluates the set of genes for DTU. +All aligned samples are analysed together with `bambu` to produce the primary +cohort transcriptome, transcript counts, gene counts, and the `RDS` objects used +for downstream differential analysis. This shared model is the main cohort-level +result and is published under `cohort/`. + +### 6. Independent per-sample transcriptomes + +Each sample is also processed separately with `bambu` so the workflow produces +sample-specific GTF, FASTA, count tables, and metadata under +`samples//`. These per-sample outputs are useful for inspecting sample +specific transcript models without changing the shared cohort transcriptome used +for DE/DTU. + +### 7. Transcript sequence generation and QC + +Transcript FASTA files are derived from GTF plus genome using `gffread`. +When `--skip_sqanti` is not set, `SQANTI3` classifies the cohort and per-sample +transcriptomes and produces structural QC summaries. The cohort `SQANTI3` +results live under `cohort/sqanti_cohort/`, while per-sample `SQANTI3` +directories are published under `samples//_sqanti/`. + +### 8. Optional DE and DTU analysis + +When `--de_analysis` is enabled, the workflow checks the experimental design, +runs `DESeq2` for differential gene expression, and runs `DEXSeq` for +differential transcript usage. These analyses use the shared `bambu` outputs +and the design columns in the sample sheet, and each comparison is written to +its own subdirectory under `de_analysis//`. + +### 9. What you need to provide + +The workflow's analysis is controlled by a user provided genome, annotation, and +`bambu` mode. + +* use `--transcriptome_mode` to choose between `discover` and + `fixed_annotation` +* `--transcriptome_source` has been removed; use `--transcriptome_mode` instead +* both `--ref_genome` and `--ref_annotation` are required in both modes +* `--ref_transcriptome` has been removed; if you want annotation-based + quantification, use `--transcriptome_mode fixed_annotation` together with + `--ref_genome` and `--ref_annotation` +* `--cdna_preprocess` enables the `pychopper` cDNA preprocessing stage, and + `--cdna_kit`, `--pychopper_backend`, and `--pychopper_opts` control that stage +* when `--de_analysis` is enabled, the sample sheet must contain `alias`, the + primary condition column, and any requested columns named in `--covariates` +* `--cdna_preprocess` must not be combined with `--direct_rna` + +### 10. How to read the output folder + +The published outputs are organised around a small number of top-level +directories: + ++ `ingress_results//` contains prepared reads, read statistics, sample + metadata, and optional `pychopper` outputs for each sample ++ `cohort/` contains the primary joint `bambu` transcriptome, count tables, + alignments, and optional cohort `SQANTI3` outputs ++ `samples//` contains the independent per-sample `bambu` outputs and + optional per-sample `SQANTI3` outputs ++ `de_analysis//` contains DE and DTU results for each contrast when + differential analysis is enabled ++ `igv_reference/` contains the published reference indexes used for IGV @@ -206,65 +300,63 @@ The final component of this isoform analysis is a stage-wise statistical test us | Nextflow parameter name | Type | Description | Help | Default | |--------------------------|------|-------------|------|---------| -| fastq | string | FASTQ files to use in the analysis. | This accepts one of three cases: (i) the path to a single FASTQ file; (ii) the path to a top-level directory containing FASTQ files; (iii) the path to a directory containing one level of sub-directories which in turn contain FASTQ files. In the first and second case, a sample name can be supplied with `--sample`. In the last case, the data is assumed to be multiplexed with the names of the sub-directories as barcodes. In this case, a sample sheet can be provided with `--sample_sheet`. | | -| bam | string | BAM or unaligned BAM (uBAM) files to use in the analysis. | This accepts one of three cases: (i) the path to a single BAM file; (ii) the path to a top-level directory containing BAM files; (iii) the path to a directory containing one level of sub-directories which in turn contain BAM files. In the first and second case, a sample name can be supplied with `--sample`. In the last case, the data is assumed to be multiplexed with the names of the sub-directories as barcodes. In this case, a sample sheet can be provided with `--sample_sheet`. | | -| transcriptome_source | string | Select how the transcriptome used for analysis should be prepared. | For differential expression analysis, use of an existing transcriptome may be preferred and so 'precomputed' should be selected. In this case the 'ref_transcriptome' parameter should be specified. To create a reference transcriptome using an existing reference genome, select 'reference guided' and specify the 'ref_genome' parameter. | reference-guided | -| ref_genome | string | Path to reference genome sequence [.fa/.fq/.fa.gz/fq.gz]. Required for reference-based workflow. | A reference genome is required for reference-based assembly of a transcriptome. | | -| ref_transcriptome | string | Transcriptome reference file. Required for precomputed transcriptome calculation and for differential expression analysis. | A reference transcriptome related to the sample under study. Must be supplied when the 'Transcriptome source' parameter has been set to 'precomputed' or to perform differential expression. | | -| ref_annotation | string | A reference annotation in GFF2 or GFF3 format (extensions .gtf(.gz), .gff(.gz), .gff3(.gz)). Only annotation files from [Encode](https://www.encodeproject.org), [Ensembl](https://www.ensembl.org/index.html) and [NCBI](https://www.ncbi.nlm.nih.gov/) are supported. | This will be used for guiding the transcriptome assembly and to label transcripts with their corresponding gene identifiers. Note: If in de_analysis mode transcript strands must be only + or -. | | -| direct_rna | boolean | Set to true for direct RNA sequencing. | Omits the pychopper step. | False | -| analyse_unclassified | boolean | Analyse unclassified reads from input directory. By default the workflow will not process reads in the unclassified directory. | If selected and if the input is a multiplex directory the workflow will also process the unclassified directory. | False | +| fastq | string | FASTQ reads to analyse. | You can provide a single FASTQ, a folder of FASTQs, or a multiplexed folder containing one sub-folder per sample or barcode. | | +| bam | string | BAM or uBAM reads to analyse. | You can provide a single BAM or uBAM, a folder of BAMs, or a multiplexed folder containing one sub-folder per sample or barcode. | | +| analyse_unclassified | boolean | Include unclassified reads from multiplexed input directories. | | False | +| analyse_fail | boolean | Include fail reads from multiplexed input directories. | | False | +| fastq_chunk | integer | Maximum number of reads per ingress chunk. | Useful mainly for testing or for splitting very large inputs into smaller pieces. | | -### Output Options +### Reference Options | Nextflow parameter name | Type | Description | Help | Default | |--------------------------|------|-------------|------|---------| -| out_dir | string | Directory for output of all user-facing files. | | output | -| igv | boolean | Visualize outputs in the EPI2ME IGV visualizer. | Enabling this option will visualize the output alignment files in the EPI2ME Desktop App IGV visualizer. | False | +| ref_genome | string | Reference genome FASTA. | Required in both discover and fixed_annotation modes. | | +| ref_annotation | string | Reference transcript annotation in GTF or GFF format. | Required in both discover and fixed_annotation modes. | | +| transcriptome_mode | string | How bambu should prepare the transcriptome model. | Use discover for reference-guided transcript discovery and quantification, or fixed_annotation for quantification only against the supplied annotation. | discover | +| direct_rna | boolean | Set this for direct RNA sequencing libraries. | | False | +| cdna_preprocess | boolean | Apply the optional cDNA preprocessing stage before alignment. | This is only relevant for cDNA libraries and must not be combined with `--direct_rna`. | False | ### Sample Options | Nextflow parameter name | Type | Description | Help | Default | |--------------------------|------|-------------|------|---------| -| sample_sheet | string | A CSV file used to map barcodes to sample aliases. The sample sheet can be provided when the input data is a directory containing sub-directories with FASTQ files. If you are running the differential expression workflow, there must be an additional column `condition` with two labels, one of which must be `control` (e.g. `control` and `treated`). Control will indicate which samples will be used as the reference. There should be at least 3 repeats for each condition. | The sample sheet is a CSV file with, minimally, columns named `barcode` and `alias`. Extra columns are allowed. | | -| sample | string | A single sample name for non-multiplexed data. Permissible if passing a single .fastq(.gz) file or directory of .fastq(.gz) files. | | | +| sample_sheet | string | CSV file describing barcodes, aliases, and optional experimental design columns. | For multiplexed runs, the sample sheet should contain both barcode and alias. For differential analysis it must also contain alias, the condition column, and any extra columns named in `--covariates`. | | +| sample | string | Single sample name for singleplexed input or to restrict multiplexed analysis to one sample. | | | -### Options for reference-based workflow +### Analysis Options | Nextflow parameter name | Type | Description | Help | Default | |--------------------------|------|-------------|------|---------| -| plot_gffcmp_stats | boolean | Create a PDF of plots from showing gffcompare results | If set to true, a PDF file containing detailed gffcompare reults will be output | True | -| gffcompare_opts | string | Extra command-line options to give to gffcompare -r | For a list of possible options see [gffcompare](https://ccb.jhu.edu/software/stringtie/gffcompare.shtml). | -R | -| minimap2_index_opts | string | Extra command-line options for minimap2 indexing. | See [minimap2 index options](https://lh3.github.io/minimap2/minimap2.html#4) for more information. These will only be relevant in the reference based transcriptome assembly. | -k 14 | -| minimap2_opts | string | Additional command-line options for minimap2 alignment. | See [minimap2 options](https://lh3.github.io/minimap2/minimap2.html#5) for further information. These will only be relevant in the reference based transcriptome assembly. | -uf | -| minimum_mapping_quality | integer | filter aligned reads by MAPQ quality. | Reads that do not meet this mapping quality after minimap2 alignment, will be filtered out. | 40 | -| stringtie_opts | string | Extra command-line options for stringtie transcript assembly. | For additional String tie options see [here](https://github.com/gpertea/stringtie#stringtie-options). | --conservative | +| de_analysis | boolean | Run differential gene expression and differential transcript usage analyses. | | False | +| condition_column | string | Main comparison column in the sample sheet. | | condition | +| covariates | string | Comma-separated extra sample-sheet columns to adjust for, for example batch. | Each listed name must exist as a column in the sample sheet. | | +| reference_level | string | Baseline group for the main comparison column. | If omitted, the workflow will use control when that level exists. | | -### Differential Expression Options +### Output Options | Nextflow parameter name | Type | Description | Help | Default | |--------------------------|------|-------------|------|---------| -| de_analysis | boolean | Run DE anaylsis | Running this requires you to provide at least two replicates for a control and treated sample as well as a sample sheet param. | False | -| min_gene_expr | integer | The minimum number of total mapped sequence reads required for a gene to be considered in differential transcript usage analysis. | Filtering at the gene level ensures that the observed transcript ratios are calculated with a minimum number of counts per gene. | 10 | -| min_feature_expr | integer | The minimum number of reads assigned to a transcript for it to be considered in differential transcript usage analysis. | Filter out transcripts that do not have this minimum number of transcript expression, reducing noise. | 3 | -| min_samps_gene_expr | integer | Set the minimum number of samples in which a gene is expressed to be included in the differential transcript usage analysis. | A gene must be expressed in at least this number of samples for the gene be included in the differential transcript usage analysis. Filtering at the gene level improves the reliability of the observed transcript ratios. | 3 | -| min_samps_feature_expr | integer | Set the minimum number of samples in which a transcript is expressed to be included in the differential transcript usage analysis. | A transcript must expressed in at least this minimum number of samples to be included in the analysis. Should be equal to the number of replicates per sample you have. | 1 | +| out_dir | string | Directory for user-facing workflow outputs. | | output | +| igv | boolean | Generate an IGV configuration file for the aligned BAM outputs. | | False | ### Advanced Options | Nextflow parameter name | Type | Description | Help | Default | |--------------------------|------|-------------|------|---------| -| threads | integer | Number of CPU threads. | Only provided to processes including alignment and and assembly that benefit from multiple threads. | 4 | -| cdna_kit | string | If cDNA reads are used, select the kit used. | This will be used by pychopper to preprocess the reads for downstream analysis. | SQK-PCS109 | -| pychopper_backend | string | Pychopper can use one of two available backends for identifying primers in the raw reads | 'edlib' is set by default due to its high performance. However, it may be less sensitive than 'phmm'. | edlib | -| pychopper_opts | string | Extra pychopper opts | See available options (here)[https://github.com/epi2me-labs/pychopper#usage] | | -| bundle_min_reads | integer | Minimum size of bam bundle for parallel processing. | | 50000 | -| isoform_table_nrows | integer | Maximum rows to dispay in the isoform report table | | 5000 | +| threads | integer | Thread count to use for the core workflow processes. | | 4 | +| minimap2_opts | string | Extra command-line options to pass to minimap2. | | | +| ndr | number | Optional bambu novel discovery rate override. | | | +| cdna_kit | string | ONT cDNA kit identifier used for pychopper preprocessing. | Only used when `--cdna_preprocess` is enabled. The workflow derives the pychopper kit code from this value. | SQK-PCS109 | +| pychopper_backend | string | Primer-detection backend to use for pychopper preprocessing. | | edlib | +| pychopper_opts | string | Extra command-line options to pass to pychopper. | | | +| skip_sqanti | boolean | Skip SQANTI3 transcript classification and QC. | | False | +| sqanti_skip_orf | boolean | Skip ORF prediction during SQANTI3 QC. | | True | +| sqanti_extra_args | string | Extra command-line options to pass to SQANTI3. | | | @@ -277,35 +369,35 @@ Output files may be aggregated including information for all samples or provided | Title | File path | Description | Per sample or aggregated | |-------|-----------|-------------|--------------------------| -| workflow report | wf-transcriptomes-report.html | a HTML report document detailing the primary findings of the workflow | aggregated | -| Per file read stats | fastq_ingress_results/{{ alias }}//reads/fastcat_stats/per-file-stats.tsv | A TSV with per file read stats, including all samples. | aggregated | -| Read stats | fastq_ingress_results/{{ alias }}//reads/fastcat_stats/per-read-stats.tsv | A TSV with per read stats, including all samples. | aggregated | -| Run ID's | fastq_ingress_results/{{ alias }}//reads/fastcat_stats/run_ids | List of run IDs present in reads. | aggregated | -| Meta map json | fastq_ingress_results/{{ alias }}//reads/metamap.json | Metadata used in workflow presented in a JSON. | aggregated | -| Concatenated sequence data | fastq_ingress_results/{{ alias }}//reads/{{ alias }}.fastq.gz | Per sample reads concatenated in to one FASTQ file. | per-sample | -| Assembled transcriptome | {{ alias }}_transcriptome.fas | Per sample assembled transcriptome. Not output if a reference annotation was supplied | per-sample | -| Annotated assembled transcriptome | {{ alias }}_merged_transcriptome.fas | Per sample annotated assembled transcriptome. Only output if a reference annotation was supplied | per-sample | -| Alignment summary statistics | {{ alias }}_read_aln_stats.tsv | Per sample alignment summary statistics. | per-sample | -| GFF compare results. | {{ alias }}_gffcompare | All GFF compare output files. | per-sample | -| Differential gene expression results | de_analysis/results_dge.tsv | This is a gene-level result file that describes genes and their probability of showing differential expression between experimental conditions. | aggregated | -| Differential gene expression report | de_analysis/results_dge.pdf | Summary report of differential gene expression analysis as a PDF. | aggregated | -| Differential transcript usage gene TSV | de_analysis/results_dtu_gene.tsv | This is a gene-level result file from DEXSeq that lists annotated genes and their probabilities of differential expression. | aggregated | -| Differential transcript usage report | de_analysis/results_dtu.pdf | Summary report of differential transcript usage results as a PDF. | aggregated | -| Differential transcript usage TSV | de_analysis/results_dtu_transcript.tsv | This is a transcript-level result file from DEXSeq that lists annotated genes and their probabilities of differential expression. | aggregated | -| Differential transcript usage stageR TSV | de_analysis/results_dtu_stageR.tsv | This is the output from StageR and it shows both gene and transcript probabilities of differential expression | aggregated | -| Differential transcript usage DEXSeq TSV | de_analysis/results_dexseq.tsv | The complete output from the DEXSeq-analysis, shows both gene and transcript probabilities of differential expression. | aggregated | -| Gene counts | de_analysis/all_gene_counts.tsv | Raw gene counts created by the Salmon tool, before filtering. | aggregated | -| Gene counts per million | de_analysis/cpm_gene_counts.tsv | This file shows counts per million (CPM) of the raw gene counts to facilitate comparisons across samples. | aggregated | -| Transcript counts | de_analysis/unfiltered_transcript_counts_with_genes.tsv | Raw transcript counts created by the Salmon tool, before filtering. Includes reference to the associated gene ID. | aggregated | -| Transcript per million counts | de_analysis/unfiltered_tpm_transcript_counts.tsv | This file shows transcripts per million (TPM) of the raw counts to facilitate comparisons across samples. | aggregated | -| Transcript counts filtered | de_analysis/filtered_transcript_counts_with_genes.tsv | Filtered transcript counts, used for differential transcript usage analysis. Includes a reference to the associated gene ID. | aggregated | -| Transcript info table | {{ alias }}_transcripts_table.tsv | This file details each isoform that was reconstructed from the input reads. It contains a subset of columns from the .tmap output from [gffcompare](https://ccb.jhu.edu/software/stringtie/gffcompare.shtml) | per-sample | -| Final non redundant transcriptome | de_analysis/final_non_redundant_transcriptome.fasta | Transcripts that were used for differential expression analysis including novel transcripts with the identifiers used for DE analysis. Only applicable when the ref_transcriptome parameter is not provided. | aggregated | -| Index of reference FASTA file | igv_reference/{{ ref_genome_file }}.fai | Reference genome index of the FASTA file required for IGV config. | aggregated | -| GZI index of the reference FASTA file | igv_reference/{{ ref_genome_file }}.gzi | GZI Index of the reference FASTA file. | aggregated | -| JSON configuration file for IGV browser | igv.json | JSON configuration file to be loaded in IGV for visualising alignments against the reference. | aggregated | -| BAM file (minimap2) | BAMS/{{ alias }}.reads_aln_sorted.bam | BAM file generated from mapping input reads to the reference. | per-sample | -| BAM index file (minimap2) | BAMS/{{ alias }}.reads_aln_sort.bam.bai | Index file generated from mapping input reads to the reference. | per-sample | +| Workflow report | wf-transcriptomes-report.html | HTML report summarising transcript discovery, quantification, optional SQANTI3 classification, and optional differential analysis results. | aggregated | +| Per-file read stats | ingress_results/{{ alias }}/fastcat_stats/per-file-stats.tsv | Read statistics for each input FASTQ file in a sample, when FASTQ read stats are available. | per-sample | +| Per-read stats | ingress_results/{{ alias }}/fastcat_stats/per-read-stats.tsv.gz | Read statistics for individual reads in a sample, when this output is enabled. | per-sample | +| Ingress reads | ingress_results/{{ alias }}/seqs.fastq.gz | Reads prepared from the input data for downstream analysis. | per-sample | +| Ingress metadata | ingress_results/{{ alias }}/metamap.json | Per-sample metadata used by the workflow. | per-sample | +| 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 | +| 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 | +| Cohort gene counts | cohort/gene_counts.tsv | Gene-level count matrix derived from bambu output. | aggregated | +| Cohort transcript metadata | cohort/transcript_metadata.tsv | Transcript annotations and bambu transcript classes for the cohort model. | aggregated | +| Cohort SQANTI3 summary | cohort/sqanti_cohort/classification_summary.tsv | SQANTI3 classification summary for the cohort transcriptome when SQANTI3 QC is enabled. | aggregated | +| Per-sample transcriptome GTF | samples/{{ alias }}/transcripts.gtf | Independent bambu transcript model for an individual sample. | per-sample | +| Per-sample transcriptome FASTA | samples/{{ alias }}/{{ alias }}.transcriptome.fa | Transcript sequences derived from the per-sample GTF. | per-sample | +| Per-sample transcript counts | samples/{{ alias }}/transcript_counts.tsv | Transcript-level abundance estimates for the per-sample bambu model. | per-sample | +| Per-sample gene counts | samples/{{ alias }}/gene_counts.tsv | Gene-level abundance estimates for the per-sample bambu model. | per-sample | +| Per-sample transcript metadata | samples/{{ alias }}/transcript_metadata.tsv | Transcript annotations and bambu transcript classes for the per-sample model. | per-sample | +| Per-sample SQANTI3 summary | samples/{{ alias }}/{{ alias }}_sqanti/classification_summary.tsv | SQANTI3 classification summary for the per-sample transcriptome when SQANTI3 QC is enabled. | per-sample | +| Differential gene expression results | de_analysis/{{ contrast }}/results_dge.tsv | DESeq2 gene-level differential expression results for one contrast. | aggregated | +| Differential gene expression plots | de_analysis/{{ contrast }}/results_dge.pdf | PDF plots generated during DESeq2 analysis for one contrast. | aggregated | +| Differential transcript usage results | de_analysis/{{ contrast }}/results_dtu_transcript.tsv | Transcript-level DTU results for one contrast. | aggregated | +| Differential transcript usage gene summary | de_analysis/{{ contrast }}/results_dtu_gene.tsv | Gene-level DTU summary for one contrast. | aggregated | +| DEXSeq results | de_analysis/{{ contrast }}/results_dexseq.tsv | Full DEXSeq result table for one contrast. | aggregated | +| Differential transcript usage plots | de_analysis/{{ contrast }}/results_dtu.pdf | PDF plots generated during DEXSeq analysis for one contrast. | aggregated | +| IGV configuration | igv.json | JSON configuration for viewing the aligned BAMs in IGV. | aggregated | +| Reference FASTA index | igv_reference/{{ ref_genome_file }}.fai | FAI index for the reference genome published for IGV. | aggregated | +| Reference GZI index | igv_reference/{{ ref_genome_file }}.gzi | GZI index for a compressed reference genome published for IGV. | aggregated | @@ -322,27 +414,216 @@ Find related RNA and cDNA sequencing protocols in the ## Troubleshooting -+ If the workflow fails please run it with the demo data set to ensure the workflow itself is working. This will help us determine if the issue is related to the environment, input parameters or a bug. -+ See how to interpret some common nextflow exit codes [here](https://labs.epi2me.io/trouble-shooting/). -+ Renaming, moving or deleting the input BAM, reference genome or the output directory from the location provided at runtime will stop IGV in the EPI2ME Desktop app from loading. ++ Check that the reference genome and reference annotation use overlapping + sequence names. The workflow checks this early and will fail if there is no + overlap at all. ++ If `--de_analysis` is enabled, ensure the sample sheet contains `alias`, the + primary condition column, and any columns named in `--covariates`. ++ DE/DTU requires at least two condition levels and at least two samples per + level. ++ `--cdna_preprocess` must not be combined with `--direct_rna`. ++ See how to interpret common Nextflow exit codes + [here](https://labs.epi2me.io/trouble-shooting/). + +### Common confusion when coming from the previous workflow version + +#### I supplied `--ref_transcriptome`, but the workflow still built or used `bambu` outputs + +The previous workflow version used `--ref_transcriptome` as a main driver for +transcript-level downstream analysis. In the current version the +`--ref_transcriptome` option has been removed and is not part of the current +`bambu` input setup. + +Use `--transcriptome_mode fixed_annotation` together with +`--ref_genome` and `--ref_annotation` if you want annotation-driven quantification +without transcript discovery. + +#### I omitted `--ref_annotation` because I expected the old input setup + +The previous workflow version could be driven from a different combination of +transcriptome inputs. In the current version both `--ref_genome` and +`--ref_annotation` are required in `discover` and `fixed_annotation` modes. +Always provide a compatible genome FASTA and transcript annotation when +launching the workflow. + +#### I used `--transcriptome_source` and got behaviour I did not expect + +In the previous workflow version the `--transcriptome_source` was the main +setting that chose how the workflow behaved. In the current version, +`--transcriptome_mode` is the main setting that controls this. The +`--transcriptome_source` option has been removed from the workflow interface. + +The options `--transcriptome_mode discover` or `--transcriptome_mode +fixed_annotation` should be used to choose between the modes of operation. + +#### I expected pychopper-style preprocessing for cDNA libraries + +cDNA preprocessing is available through +`--cdna_preprocess`, which runs `pychopper` before alignment. The +parameters for that stage are `--cdna_kit`, `--pychopper_backend`, and +`--pychopper_opts`. + +Enable `--cdna_preprocess` for cDNA libraries when you want +`pychopper` preprocessing, or use `--direct_rna` without +`--cdna_preprocess` for direct RNA data. + +#### I cannot find the old flat DE output files + +In the previous workflow version, DE files appeared directly under `de_analysis/`. +In the current version DE and DTU results are grouped by contrast under +`de_analysis//`. Look for outputs such as +`de_analysis//results_dge.tsv` and +`de_analysis//results_dtu_transcript.tsv`. + +#### I expected the old output layout or transcriptome files + +The previous workflow version emitted one flat set of transcriptome. +In the current versions, the output folder is organised around `ingress_results/`, +`cohort/`, `samples//`, `de_analysis//`, and +`igv_reference/`. + +Primary shared transcriptome results are under `cohort/`, and +`samples//` for sample-specific models. If `--cdna_preprocess` is enabled, +look under `ingress_results//` for the `pychopper` outputs associated +with each sample. + +#### My DE/DTU run fails because of `--sample_sheet` + +The sample sheet rules have been amended compared to the previous version +to allow for multi-way comparisons. + +The sample sheet must contain `alias`, the primary condition +column, and any columns named in `--covariates`. At least two condition levels +are required, and each level must contain at least two samples. + +What to change: verify the sample sheet columns first, then check +`--condition_column`, `--covariates`, and `--reference_level`. + +#### I combined `--direct_rna` with `--cdna_preprocess` + +This is not a valid combination of parameters. Use `--direct_rna` for direct RNA +libraries and leave `--cdna_preprocess` unset. + +#### I hit genome/annotation validation or strand-related annotation warnings + +The workflow validates that the annotation and genome share +sequence names, and it excludes unstranded annotation entries from the +differential analysis path. + +Confirm the genome and annotation come from a compatible source, +and ensure the annotation uses only `+` or `-` strand values where required for +DE/DTU and `SQANTI3`. + ## FAQs -*Does the workflow support de novo assembly?* - Currently the workflow does not have a *de novo* mode. +### Does the workflow support both cDNA and direct RNA? -*Why is the IGV panel not showing?* - The workflow expects either an uncompressed or [`bgzip`](https://www.htslib.org/doc/bgzip.html)-compressed reference. If the user provides a reference compressed not with `bgzip`, the workflow will run to completion, but won't be able to generate the necessary indexes to visualize the outputs in IGV. +Yes. Use `--direct_rna` for direct RNA data. cDNA is the default mode. + +### Do I need both `--ref_genome` and `--ref_annotation` in fixed-annotation mode? + +Yes. `bambu` still uses the genome together with the imported annotation. + +### Does the workflow create one shared transcriptome or one per sample? + +Both. The joint cohort model is the primary result for reporting and DE/DTU, and +the workflow also emits independent per-sample transcriptomes. + +### Can I run DE/DTU without transcript discovery? + +Yes. Use `--transcriptome_mode fixed_annotation` together with `--de_analysis`. + +### What changed from the previous workflow version? + +The current workflow uses `bambu`, optional `SQANTI3`, `DESeq2`, and `DEXSeq`. +The main transcriptome result is now one shared `bambu` model built from all +samples together, with separate per-sample `bambu` outputs published alongside +it. The most important differences are summarised below. + +#### Why do the results now mention `bambu` and `SQANTI3` instead of StringTie or GffCompare? + +The workflow now uses a different set of transcript analysis tools. It builds +its transcript models with `bambu`, optionally classifies them with `SQANTI3`, +and performs DE/DTU from the cohort `bambu` outputs. + +#### Why is `--ref_annotation` now required even in fixed-annotation mode? + +`bambu` still requires the annotation together with the genome in both +`discover` and `fixed_annotation` modes. Fixed-annotation mode means +annotation-driven quantification, not annotation-only execution. + +#### Can I still use `--ref_transcriptome`? + +No. `--ref_transcriptome` has been removed from the workflow interface. + +Short old-to-new example: + +```text +Previous workflow version: + --transcriptome_source precomputed --ref_transcriptome transcripts.fa + +Current workflow: + --transcriptome_mode fixed_annotation \ + --ref_genome genome.fa \ + --ref_annotation annotation.gtf +``` + +#### Why do I now get both cohort and per-sample transcriptomes? + +The workflow now treats both as important outputs. The shared cohort model is +the main transcriptome used for reporting and optional DE/DTU, while the +per-sample transcriptomes are provided for looking at each sample separately. + +#### Why are DE results under `de_analysis//`? + +The workflow now writes one subdirectory per comparison instead of publishing +one flat DE result set. This makes the output folder clearer when more than one +comparison is present. + +Short old-to-new example: + +```text +Previous workflow version: + de_analysis/results_dge.tsv + +Current workflow: + de_analysis//results_dge.tsv +``` + +#### Can I still run fixed-annotation quantification without transcript discovery? + +Yes. Use `--transcriptome_mode fixed_annotation` together with `--ref_genome`, +`--ref_annotation`, and any optional DE/DTU settings. + +#### What should I expect to differ in report contents and output filenames? + +Expect the report and output folder to emphasise: + ++ optional `pychopper` preprocessing outputs under `ingress_results//` ++ the joint cohort `bambu` transcriptome under `cohort/` ++ the per-sample `bambu` transcriptomes under `samples//` ++ optional `SQANTI3` results under cohort and per-sample directories ++ contrast-specific DE/DTU outputs under `de_analysis//` + +If your question is not answered here, please report issues or suggestions on +the [GitHub issues](https://github.com/epi2me-labs/wf-transcriptomes/issues) +page or start a discussion on the +[community](https://community.nanoporetech.com/). -If your question is not answered here, please report any issues or suggestions on the [github issues](https://github.com/epi2me-labs/wf-transcriptomes/issues) page or start a discussion on the [community](https://community.nanoporetech.com/). ## Related blog posts -+ [How to align your data](https://labs.epi2me.io/how-to-align/) ++ See the [EPI2ME website](https://epi2me.nanoporetech.com/) for more workflow resources and blog posts. -See the [EPI2ME website](https://labs.epi2me.io/) for lots of other resources and blog posts. +### References + +[Systematic assessment of long-read RNA-seq methods for transcript identification and quantification +](https://www.nature.com/articles/s41592-024-02298-3) diff --git a/bin/de_analysis.R b/bin/de_analysis.R deleted file mode 100755 index cedde4f..0000000 --- a/bin/de_analysis.R +++ /dev/null @@ -1,226 +0,0 @@ -#!/usr/bin/env Rscript - -suppressMessages(library(argparser)) - -parser <- arg_parser("Run differential expression analysis") -parser <- add_argument(parser, "--annotation", help="Reference annotation.") -parser <- add_argument(parser, "--min_samps_gene_expr", help="Minimum number of samples a gene must be expressed in to be included in differential gene expression.", type="numeric") -parser <- add_argument(parser, "--min_samps_feature_expr", help="Minimum number of samples for differential transcript usage.", type="numeric") -parser <- add_argument(parser, "--min_gene_expr", help="Minimum counts per gene required for differential gene expression.", type="numeric") -parser <- add_argument(parser, "--min_feature_expr", help="Minimum counts per transcript required for differential transcript usage.", type="numeric") -parser <- add_argument(parser, "--sample_sheet", help="Sample sheet.") -parser <- add_argument(parser, "--all_counts", help="All transcript counts CSV file.") -parser <- add_argument(parser, "--de_out_dir", help="Directory where differential expression out files will be saved. Directory will be created if it does not exist", default="de_analysis") -parser <- add_argument(parser, "--merged_out_dir", help="Directory where merged count files will be saved. Directory will be created if it does not exist", default="merged") -argv <- parse_args(parser) - -suppressMessages(library("DRIMSeq")) -suppressMessages(library("GenomicFeatures")) -suppressMessages(library("edgeR")) - -# Some functions, including dmDSdata, converts '.' in sample IDs to '-'. -# Make the output DF match the sample IDs in the sample sheet -# start_col is the index of the first sample column in the data frame -rename_sample_columns <- function(df, alias, start_col) { - colnames(df)[start_col:ncol(df)] <- alias - return(df) -} - -# Create output directories -if (!dir.exists(argv$de_out_dir)){ - dir.create(argv$de_out_dir, recursive=TRUE) -} -if (!dir.exists(argv$merged_out_dir)){ - dir.create(argv$merged_out_dir, recursive=TRUE) -} - -cat("Loading counts, conditions and parameters.\n") -cts <- as.matrix(read.csv(argv$all_counts, sep="\t", row.names="Reference", stringsAsFactors=FALSE, check.names=FALSE)) -# Set up sample data frame: -coldata <- read.csv(argv$sample_sheet, row.names="alias", sep=",", stringsAsFactors=TRUE, check.names=FALSE) -coldata$alias <- rownames(coldata) -#dmDSdata looks up sample_id in coldata, we need this to be sample alias -coldata$sample_id <- coldata$alias - -# Reorder the input counts columns to match the sample sheet order. -# This ensures we don't missassign column names, when renaming output DF columns. -cts <- cts[, coldata$alias, drop = FALSE] -# check if control condition exists, sets as reference -if(!"control" %in% coldata$condition) - stop("sample_sheet.csv does not contain 'control' - condition - unable to set reference.") -coldata$condition <- relevel(coldata$condition, ref = "control") - -# a .gff annotation file extension may be gff2(gtf) or gff3 so check in files for use of = in the attribute field -# if '=' present it is gff3 if not it is gtf. -# see https://www.ensembl.org/info/website/upload/gff.html -# and http://gmod.org/wiki/GFF2#Converting_GFF2_to_GFF3 -cat("Checking annotation file type.\n") -lines <- readLines(file(argv$annotation), n=10000) -# If transcript_id containing '=' (format eg. transcript_id=xxx) -# annotation type is gff3 -check_file_type <- sum(grepl("transcript_id=", lines)) -if (check_file_type != 0){ - cat("Annotation file type is gff3.\n") - annotation_type <- "gff3" -} else { - # otherwise gtf - cat("Annotation file type is gtf.\n") - annotation_type <- "gtf" -} - -# Transcript_id versions (eg. ENTXXX.1, eg. ENTXXX.2) represent how many times that transcript reference has been changed -# during its time in the database. -# Not all annotation files include it as part of the transcript_id - notably Ensembl -# The following handles this. -cat("Checking annotation file for presence of transcript_id versions.\n") -# Get the first transcript_id from the annotation file by parsing -lines <- readLines(file(argv$annotation), n=100000) -# Find transcript_ids in first 1000 lines and check if they contain dot (format eg. ENTXXX.1) -check_version <- sum(grepl("transcript_id[^;]+\\.", lines)) -if (check_version != 0){ - # we do not need to strip the count file rows if ref_annotation includes versions - cat("Annotation file transcript_ids include versions.\n") - } else { - # otherwise remove the versions - rownames(cts) <- lapply(rownames(cts), sub, pattern = "\\.\\d+$", replacement = "") - cat("Annotation file transcript_ids do not include versions so also strip versions from the counts df.\n") - } - -cat("Loading annotation database.\n") -txdb <- makeTxDbFromGFF(argv$annotation, format = annotation_type) -txdf <- select(txdb, keys(txdb,"GENEID"), "TXNAME", "GENEID") -tab <- table(txdf$GENEID) -txdf$ntx<- tab[match(txdf$GENEID, names(tab))] - - -cts <- cts[rownames(cts) %in% txdf$TXNAME, ] # FIXME: filter for transcripts which are in the annotation. Why they are not all there? -# Reorder transcript/gene database to match input counts: -txdf <- txdf[match(rownames(cts), txdf$TXNAME), ] -rownames(txdf) <- NULL - -# Create counts data frame: -counts<-data.frame(gene_id=txdf$GENEID, feature_id=txdf$TXNAME, cts) - -counts <- rename_sample_columns(counts, coldata$alias, 3) -# output unfiltered version of the counts table now we have paired transcripts with gene ids -write.table(counts, file=file.path(argv$de_out_dir, "unfiltered_transcript_counts_with_genes.tsv"), sep="\t", row.names = FALSE, quote=FALSE) - - -cat("Filtering counts using DRIMSeq.\n") - -d <- dmDSdata(counts=counts, samples=coldata) -trs_cts_unfiltered <- counts(d) - -d <- dmFilter(d, min_samps_gene_expr=argv$min_samps_gene_expr, min_samps_feature_expr=argv$min_samps_feature_expr, - min_gene_expr=argv$min_gene_expr, min_feature_expr=argv$min_feature_expr) - -cat("Building model matrix.\n") -design <- model.matrix(~condition, data=DRIMSeq::samples(d)) - - -suppressMessages(library("dplyr")) - -# Sum transcript counts into gene counts: -cat("Sum transcript counts into gene counts.\n") -trs_cts <- counts(d) -trs_cts <- rename_sample_columns(trs_cts, coldata$alias, 3) -write.table(trs_cts, file=file.path(argv$merged_out_dir, "filtered_transcript_counts_with_genes.tsv"), sep="\t", row.names = FALSE, quote=FALSE) - -gene_cts <- trs_cts_unfiltered %>% dplyr::select(c(1, 3:ncol(trs_cts))) %>% group_by(gene_id) %>% summarise_all(tibble::lst(sum)) %>% data.frame() -rownames(gene_cts) <- gene_cts$gene_id -gene_cts$gene_id <- NULL -gene_cts <- rename_sample_columns(gene_cts, coldata$alias, 1) -write.table(gene_cts, file=file.path(argv$merged_out_dir, "all_gene_counts.tsv"), sep="\t", quote=FALSE) - -# Output count per million of the gene counts using edgeR CPM -cpm_gene_counts <- cpm(gene_cts) -# Add gene_id as index column header -cpm_gene_counts <- cbind(var_name = rownames(cpm_gene_counts), cpm_gene_counts) -rownames(cpm_gene_counts) <- NULL -colnames(cpm_gene_counts)[1] <- "gene_id" -cpm_gene_counts <- rename_sample_columns(cpm_gene_counts, coldata$alias, 2) -write.table(cpm_gene_counts, file=file.path(argv$de_out_dir, "cpm_gene_counts.tsv"), sep="\t", quote=FALSE, row.names = FALSE) - -# Differential gene expression using edgeR: -cat("Running differential gene expression analysis using edgeR.\n") - -y <- DGEList(gene_cts) -y <- calcNormFactors(y) -y <- estimateDisp(y,design) -fit <- glmQLFit(y,design) -qlf <- glmQLFTest(fit) -edger_res <- topTags(qlf, n=nrow(y), sort.by="PValue")[[1]] - -pdf("de_analysis/results_dge.pdf") - -# create status vector -status <- ifelse( - qlf$PValue<0.01 & qlf$logFC>0, - 'up', - ifelse( - qlf$PValue<0.01 & qlf$logFC<=0, - 'down', - 'notsig' - ) -) -plotMD(qlf, status=status, values=c("up","down","notsig"), hl.col=c("red","blue","black")) -abline(h=c(-1,1), col="blue") -plotQLDisp(fit) - -write.table(as.data.frame(edger_res), file=file.path(argv$de_out_dir, "results_dge.tsv"), sep="\t") - -# Differential transcript usage using DEXSeq: -suppressMessages(library("DEXSeq")) -cat("Running differential transcript usage analysis using DEXSeq.\n") - -sample.data<-DRIMSeq::samples(d) -count.data <- round(as.matrix(counts(d)[,-c(1:2)])) -dxd <- DEXSeqDataSet(countData=count.data, sampleData=sample.data, design=~sample + exon + condition:exon, featureID=trs_cts$feature_id, groupID=trs_cts$gene_id) -dxd <- estimateSizeFactors(dxd) -dxd <- estimateDispersions(dxd) -dxd <- testForDEU(dxd, reducedModel=~sample + exon) -dxd <- estimateExonFoldChanges( dxd, fitExpToVar="condition") -dxr <- DEXSeqResults(dxd, independentFiltering=FALSE) - -dev.off() -pdf("de_analysis/results_dtu.pdf") -plotMA(dxr, cex=0.8, alpha=0.05) -plotDispEsts(dxd) - -qval <- perGeneQValue(dxr) -dxr.g<-data.frame(gene=names(qval), qval) -dxr.g <- dxr.g[order(dxr.g$qval),] - -dxr_out <- as.data.frame(dxr[,c("featureID", "groupID", "pvalue")]) -dxr_out <- dxr_out[order(dxr$pvalue),] - -write.table(dxr.g, file=file.path(argv$de_out_dir, "results_dtu_gene.tsv"), sep="\t") -write.table(dxr_out, file=file.path(argv$de_out_dir, "results_dtu_transcript.tsv"), sep="\t") - -# and writing out some of the DEXSeq metrics to accompany EPI2ME Labs tutorial -colnames(dxr)[grep("log2fold", colnames(dxr))] <- "log2fold" -MADTUdata <- data.frame(dxr)[order(dxr$padj),c("exonBaseMean", "log2fold", "pvalue", "padj")] -MADTUdata$exonBaseMean <- log2(MADTUdata$exonBaseMean) -colnames(MADTUdata)[which(colnames(MADTUdata)=="exonBaseMean")] <- "Log2MeanExon" -colnames(MADTUdata)[which(colnames(MADTUdata)=="log2fold")] <- "Log2FC" -write.table(MADTUdata, file=file.path(argv$de_out_dir, "results_dexseq.tsv"), sep="\t") - -# stageR analysis of DEXSeq results: -cat("stageR analysis\n") -library(stageR) - -cat("Running stageR analysis on the differential transcript usage results.\n") -pConfirmation <- matrix(dxr$pvalue, ncol=1) - -dimnames(pConfirmation) <- list(dxr$featureID, "transcript") -pScreen <- qval -tx2gene <- as.data.frame(dxr[,c("featureID", "groupID")]) - -stageRObj <- stageRTx(pScreen=pScreen, pConfirmation=pConfirmation, pScreenAdjusted=TRUE, tx2gene=tx2gene) -# note: the choice of 0.05 here means you can *only* threshold at 5% OFDR later -stageRObj <- stageWiseAdjustment(stageRObj, method="dtu", alpha=0.10) -suppressWarnings({dex.padj <- getAdjustedPValues(stageRObj, order=FALSE, onlySignificantGenes=FALSE)}) - -# dex.padj <- dex.padj[,-1] -write.table(dex.padj, file=file.path(argv$de_out_dir, "results_dtu_stageR.tsv"), sep="\t") diff --git a/bin/plot_dtu_results.R b/bin/plot_dtu_results.R deleted file mode 100755 index 16df20a..0000000 --- a/bin/plot_dtu_results.R +++ /dev/null @@ -1,65 +0,0 @@ -#!/usr/bin/env Rscript - -suppressMessages(library(argparser)) - -parser <- arg_parser("Plot results") -parser <- add_argument(parser, "--counts", help="Filtered transcript counts with genes.") -parser <- add_argument(parser, "--results_dtu", help="stageR results.") -parser <- add_argument(parser, "--sample_sheet", help="Sample sheet.") -parser <- add_argument(parser, "--pdf_out", help="PDF file name.") -argv <- parse_args(parser) - -suppressMessages(library(dplyr)) -suppressMessages(library(ggplot2)) -suppressMessages(library(tidyr)) - -# Set up sample data frame: -coldata <- read.csv(argv$sample_sheet, row.names="alias", sep=",") -coldata$condition <- factor(coldata$condition, levels=rev(levels(coldata$condition))) -coldata$type <-NULL -coldata$patient <-NULL - -# Read stageR results: -stageR <- read.csv(argv$results_dtu, sep="\t") -names(stageR) <- c("gene_id", "transcript_id", "p_gene", "p_transcript"); - -# Read filtered counts: -counts <- read.csv(argv$counts, sep="\t"); -names(counts)[2]<-"transcript_id" - -# Join counts and stageR results: -df <- counts %>% left_join(stageR, by = c("gene_id", "transcript_id")) -df <- df[order(df$p_gene),] - -scols <- setdiff(names(df),c("gene_id", "transcript_id", "p_gene", "p_transcript")) - -# Normalise counts: -for(sc in scols){ - df[sc] <- df[sc] / sum(df[sc]) -} - -# Melt data frame: -tdf <- df %>% gather(key='sample', value='norm_count',-gene_id, -transcript_id, -p_gene, -p_transcript) - -# Add sample group column: -sampleToGroup<-function(x){ - return(coldata[x,]$condition) -} - -tdf$group <- sampleToGroup(tdf$sample) - -# Filter for significant genes: -sig_level <- 0.05 -genes <- as.character(tdf[which(tdf$p_gene < sig_level),]$gene_id) -genes <- unique(genes) - -pdf(argv$pdf_out) - -for(gene in genes){ - gdf<-tdf[which(tdf$gene_id==gene),] - p_gene <- unique(gdf$p_gene) - dtu_plot <- ggplot(gdf, aes(x=transcript_id, y=norm_count)) + geom_bar(stat="identity", aes(fill=sample), position="dodge") - dtu_plot <- dtu_plot + facet_wrap(~ group) + coord_flip() - dtu_plot <- dtu_plot + ggtitle(paste(gene," : p_value=",p_gene,sep="")) - print(dtu_plot) -} diff --git a/bin/run_bambu.R b/bin/run_bambu.R new file mode 100755 index 0000000..ab428ef --- /dev/null +++ b/bin/run_bambu.R @@ -0,0 +1,449 @@ +#!/usr/bin/env Rscript + +# Set seed for reproducibility +set.seed(42) + +suppressPackageStartupMessages({ + library(argparser) + library(bambu) + library(Rsamtools) + library(SummarizedExperiment) + library(jsonlite) +}) + +parser <- arg_parser("Run bambu transcript discovery and quantification.") +parser <- add_argument(parser, "--bam_dir", help = "Directory containing BAM files.") +parser <- add_argument(parser, "--bam_path", help = "Path to a single BAM file.") +parser <- add_argument(parser, "--sample_alias", help = "Alias to use for a single BAM file.") +parser <- add_argument(parser, "--sample_sheet", help = "Optional sample sheet CSV.") +parser <- add_argument(parser, "--annotation", help = "Reference annotation GTF/GFF.") +parser <- add_argument(parser, "--genome", help = "Reference genome FASTA.") +parser <- add_argument(parser, "--transcriptome_mode", help = "discover or fixed_annotation.", default = "discover") +parser <- add_argument(parser, "--threads", help = "Number of worker threads.", type = "numeric", default = 1) +parser <- add_argument(parser, "--ndr", help = "Optional novel discovery rate.", type = "numeric") +parser <- add_argument(parser, "--out_dir", help = "Output directory.") +argv <- parse_args(parser) + +arg_missing <- function(value) { + if (is.null(value) || length(value) == 0 || all(is.na(value))) { + return(TRUE) + } + if (is.character(value)) { + return(all(!nzchar(value))) + } + FALSE +} + +required_args <- c("annotation", "genome", "out_dir") +missing_args <- required_args[vapply(required_args, function(arg_name) { + value <- argv[[arg_name]] + arg_missing(value) +}, logical(1))] +if (length(missing_args) > 0) { + stop(sprintf( + "Missing required arguments: %s", + paste(sprintf("--%s", missing_args), collapse = ", ") + )) +} + +if (arg_missing(argv$bam_dir) == arg_missing(argv$bam_path)) { + stop("Provide exactly one of --bam_dir or --bam_path.") +} + +dir.create(argv$out_dir, showWarnings = FALSE, recursive = TRUE) + +sample_df <- NULL +if (!arg_missing(argv$sample_sheet)) { + sample_df <- read.csv(argv$sample_sheet, check.names = FALSE, stringsAsFactors = FALSE) + if (!"alias" %in% names(sample_df)) { + stop("Sample sheet must contain an 'alias' column.") + } +} + +strip_alias <- function(path) { + name <- basename(path) + name <- sub("\\.aligned\\.sorted\\.bam$", "", name) + name <- tools::file_path_sans_ext(name) + name +} + +if (!arg_missing(argv$bam_dir)) { + bam_paths <- sort(list.files(argv$bam_dir, pattern = "\\.bam$", full.names = TRUE)) + if (length(bam_paths) < 1) { + stop("No BAM files were found in bam_dir.") + } + aliases <- vapply(bam_paths, strip_alias, character(1)) +} else { + bam_paths <- argv$bam_path + aliases <- if (!arg_missing(argv$sample_alias)) argv$sample_alias else strip_alias(argv$bam_path) +} + +if (!is.null(sample_df)) { + missing_aliases <- setdiff(aliases, sample_df$alias) + if (length(missing_aliases) > 0) { + stop(sprintf( + "Sample sheet is missing alias rows for BAM files: %s", + paste(missing_aliases, collapse = ", ") + )) + } + sample_df <- sample_df[match(aliases, sample_df$alias), , drop = FALSE] +} else { + sample_df <- data.frame(alias = aliases, stringsAsFactors = FALSE) +} + +annotation_obj <- prepareAnnotations(argv$annotation) +reads <- if (length(bam_paths) == 1) bam_paths else BamFileList(bam_paths, yieldSize = 1000000) + +# Handle NDR parameter with validation and documentation +default_ndr <- 0.1 +ndr_value <- default_ndr + +if (!arg_missing(argv$ndr)) { + if (argv$ndr < 0 || argv$ndr > 1) { + stop("NDR (Novel Discovery Rate) must be between 0 and 1") + } + ndr_value <- argv$ndr + message(sprintf("Using user-specified NDR = %.3f", ndr_value)) +} else { + message(sprintf("Using default NDR = %.3f", default_ndr)) +} + +if (identical(argv$transcriptome_mode, "discover")) { + message("Novel Discovery Rate (NDR) controls transcript discovery stringency:") + message(" Lower NDR (e.g., 0.05) = fewer false positive transcripts, may miss real ones") + message(" Higher NDR (e.g., 0.2) = more sensitive discovery, more false positives") + message(sprintf(" Current NDR = %.3f balances precision and recall", ndr_value)) +} + +bambu_args <- list( + reads = reads, + annotations = annotation_obj, + genome = argv$genome, + ncore = as.integer(argv$threads), + discovery = identical(argv$transcriptome_mode, "discover") +) + +if (identical(argv$transcriptome_mode, "discover")) { + bambu_args$NDR <- ndr_value +} + +message("Running bambu...") +se <- do.call(bambu, bambu_args) +message("Bambu completed successfully") +colnames(se) <- aliases + +counts_mat <- assays(se)$counts +full_length_mat <- assays(se)$fullLengthCounts + +# Collect QC statistics before filtering +qc_stats <- list() +qc_stats$total_transcripts_before_filter <- nrow(se) +qc_stats$total_genes_before_filter <- length(unique(rowData(se)$GENEID)) +qc_stats$samples <- ncol(se) + +# Filter low-count transcripts +if (is.null(full_length_mat)) { + keep_idx <- rowSums(counts_mat) > 0 +} else { + keep_idx <- rowSums(full_length_mat) > 0 +} +if (!any(keep_idx)) { + keep_idx <- rowSums(counts_mat) >= 0 +} + +qc_stats$transcripts_filtered <- sum(!keep_idx) +message(sprintf("Filtering: keeping %d / %d transcripts", sum(keep_idx), length(keep_idx))) + +se <- se[keep_idx, ] + +# Library size statistics and warnings +lib_sizes <- colSums(assays(se)$counts) +qc_stats$library_sizes <- as.list(lib_sizes) +qc_stats$min_library_size <- min(lib_sizes) +qc_stats$max_library_size <- max(lib_sizes) +qc_stats$median_library_size <- median(lib_sizes) + +if (length(lib_sizes) > 1) { + lib_size_ratio <- max(lib_sizes) / min(lib_sizes) + qc_stats$library_size_ratio <- lib_size_ratio + + if (lib_size_ratio > 3) { + warning(sprintf( + "Large library size variation detected (%.1fx difference).\n Min: %d, Max: %d reads.\n CPM normalization may not be appropriate for such variation.", + lib_size_ratio, min(lib_sizes), max(lib_sizes) + )) + qc_stats$library_size_warning <- sprintf("%.1fx variation (>3x threshold)", lib_size_ratio) + } +} + +# Per-sample detection statistics +qc_stats$transcripts_detected_per_sample <- as.list(colSums(assays(se)$counts > 0)) +qc_stats$median_transcripts_detected <- median(colSums(assays(se)$counts > 0)) + +qc_stats$total_transcripts_after_filter <- nrow(se) +qc_stats$total_genes_after_filter <- length(unique(rowData(se)$GENEID)) + +row_ranges <- rowRanges(se) +writeToGTF(row_ranges, file = file.path(argv$out_dir, "transcripts.gtf")) + +gene_se <- transcriptToGeneExpression(se) +colnames(gene_se) <- aliases + +saveRDS(se, file.path(argv$out_dir, "bambu_transcripts.rds")) +saveRDS(gene_se, file.path(argv$out_dir, "bambu_genes.rds")) +write.csv(sample_df, file.path(argv$out_dir, "samples.csv"), row.names = FALSE, quote = FALSE) + +tx_meta <- as.data.frame(rowData(se)) +if (!"TXNAME" %in% names(tx_meta)) { + tx_meta$TXNAME <- rownames(se) +} +if (!"GENEID" %in% names(tx_meta)) { + tx_meta$GENEID <- NA_character_ +} +gene_meta <- as.data.frame(rowData(gene_se)) +if (!"GENEID" %in% names(gene_meta)) { + gene_meta$GENEID <- rownames(gene_se) +} + +matrix_to_df <- function(se_obj, assay_name, id_col, meta_df) { + assay_df <- as.data.frame(assays(se_obj)[[assay_name]]) + assay_df[[id_col]] <- rownames(se_obj) + assay_df <- assay_df[, c(id_col, setdiff(names(assay_df), id_col)), drop = FALSE] + merge(meta_df, assay_df, by.x = id_col, by.y = id_col, all.y = TRUE, sort = FALSE) +} + +normalise_tsv_value <- function(value) { + if (length(value) == 0 || all(is.na(value))) { + return(NA_character_) + } + if (is.list(value)) { + value <- unlist(value, recursive = TRUE, use.names = FALSE) + } + if (length(value) == 0 || all(is.na(value))) { + return(NA_character_) + } + paste(as.character(value), collapse = ";") +} + +normalise_tsv_df <- function(df) { + as.data.frame( + lapply(df, function(column) { + if (is.list(column)) { + vapply(column, normalise_tsv_value, character(1)) + } else { + column + } + }), + stringsAsFactors = FALSE, + check.names = FALSE + ) +} + +extract_gtf_attribute <- function(attr_field, key) { + match <- regexec(sprintf('%s "([^"]*)";', key), attr_field, perl = TRUE) + captures <- regmatches(attr_field, match)[[1]] + if (length(captures) < 2) { + return(NULL) + } + captures[2] +} + +normalise_gtf_attribute_value <- function(value) { + if (is.null(value)) { + return(NULL) + } + value <- gsub('[";]', "", value) + value <- trimws(gsub("\\s+", " ", value)) + if (!nzchar(value)) { + return(NULL) + } + value +} + +sanitise_gtf_file <- function(path) { + lines <- readLines(path, warn = FALSE) + cleaned_lines <- vapply(lines, function(line) { + if (!nzchar(line) || startsWith(line, "#")) { + return(line) + } + + fields <- strsplit(line, "\t", fixed = TRUE)[[1]] + if (length(fields) < 9) { + return(line) + } + + attr_field <- fields[9] + transcript_id <- normalise_gtf_attribute_value( + extract_gtf_attribute(attr_field, "transcript_id") + ) + gene_id <- normalise_gtf_attribute_value( + extract_gtf_attribute(attr_field, "gene_id") + ) + + if (!is.null(gene_id) && grepl("\\btranscript_id\\b", gene_id)) { + gene_id <- transcript_id + } + if (is.null(gene_id)) { + gene_id <- transcript_id + } + + if (!is.null(gene_id)) { + attr_field <- sub( + 'gene_id "([^"]*)";', + sprintf('gene_id "%s";', gene_id), + attr_field, + perl = TRUE + ) + } + if (!is.null(transcript_id)) { + attr_field <- sub( + 'transcript_id "([^"]*)";', + sprintf('transcript_id "%s";', transcript_id), + attr_field, + perl = TRUE + ) + } + + fields[9] <- attr_field + paste(fields, collapse = "\t") + }, character(1)) + + writeLines(cleaned_lines, path) +} + +tx_meta <- normalise_tsv_df(tx_meta) +gene_meta <- normalise_tsv_df(gene_meta) + +write.table( + tx_meta, + file = file.path(argv$out_dir, "transcript_metadata.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE +) +write.table( + gene_meta, + file = file.path(argv$out_dir, "gene_metadata.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE +) + +tx_counts <- matrix_to_df(se, "counts", "TXNAME", tx_meta) +tx_cpm <- matrix_to_df(se, "CPM", "TXNAME", tx_meta) +gene_counts <- matrix_to_df(gene_se, "counts", "GENEID", gene_meta) +gene_cpm <- matrix_to_df(gene_se, "CPM", "GENEID", gene_meta) + +write.table( + tx_counts, + file = file.path(argv$out_dir, "transcript_counts.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE +) +write.table( + tx_cpm, + file = file.path(argv$out_dir, "transcript_cpm.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE +) +write.table( + gene_counts, + file = file.path(argv$out_dir, "gene_counts.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE +) +write.table( + gene_cpm, + file = file.path(argv$out_dir, "gene_cpm.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE +) + +sanitise_gtf_file(file.path(argv$out_dir, "transcripts.gtf")) + +# Write QC statistics as JSON for HTML report +qc_stats$transcriptome_mode <- argv$transcriptome_mode +qc_stats$ndr_used <- if (identical(argv$transcriptome_mode, "discover")) ndr_value else "N/A" +qc_stats$timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S") + +write_json( + qc_stats, + file.path(argv$out_dir, "bambu_qc_stats.json"), + pretty = TRUE, + auto_unbox = TRUE +) + +format_count <- function(value) { + if (length(value) == 0 || all(is.na(value))) { + return("NA") + } + format( + round(as.numeric(value), 0), + scientific = FALSE, + trim = TRUE, + big.mark = "," + ) +} + +# Write human-readable QC summary +qc_summary <- c( + "Bambu Quantification QC Summary", + "================================", + "", + sprintf("Timestamp: %s", qc_stats$timestamp), + sprintf("Mode: %s", argv$transcriptome_mode), + if (identical(argv$transcriptome_mode, "discover")) sprintf("NDR: %.3f", ndr_value) else NULL, + "", + "Sample Statistics:", + sprintf(" Samples analyzed: %s", format_count(qc_stats$samples)), + sprintf( + " Median library size: %s reads", + format_count(qc_stats$median_library_size) + ), + sprintf( + " Library size range: %s - %s reads", + format_count(qc_stats$min_library_size), + format_count(qc_stats$max_library_size) + ), + if (!is.null(qc_stats$library_size_warning)) sprintf(" WARNING: %s", qc_stats$library_size_warning) else NULL, + "", + "Transcript Discovery:", + sprintf( + " Transcripts before filtering: %s", + format_count(qc_stats$total_transcripts_before_filter) + ), + sprintf( + " Transcripts after filtering: %s", + format_count(qc_stats$total_transcripts_after_filter) + ), + sprintf( + " Transcripts removed: %s", + format_count(qc_stats$transcripts_filtered) + ), + sprintf( + " Median transcripts detected per sample: %s", + format_count(qc_stats$median_transcripts_detected) + ), + "", + "Gene-Level Summary:", + sprintf( + " Unique genes (before filter): %s", + format_count(qc_stats$total_genes_before_filter) + ), + sprintf( + " Unique genes (after filter): %s", + format_count(qc_stats$total_genes_after_filter) + ), + "" +) + +writeLines(qc_summary, file.path(argv$out_dir, "bambu_qc_summary.txt")) +message("QC statistics written to bambu_qc_stats.json and bambu_qc_summary.txt") + +# Save session info for reproducibility +writeLines(capture.output(sessionInfo()), file.path(argv$out_dir, "session_info.txt")) +message("Session info saved for reproducibility") diff --git a/bin/run_de_analysis.R b/bin/run_de_analysis.R new file mode 100755 index 0000000..9bdd12b --- /dev/null +++ b/bin/run_de_analysis.R @@ -0,0 +1,717 @@ +#!/usr/bin/env Rscript + +# Set seed for reproducibility +set.seed(42) + +suppressPackageStartupMessages({ + library(argparser) + library(DESeq2) + library(DEXSeq) + library(SummarizedExperiment) + library(jsonlite) +}) + +parser <- arg_parser("Run DESeq2 and DEXSeq on bambu output.") +parser <- add_argument(parser, "--transcript_rds", help = "bambu transcript RDS.") +parser <- add_argument(parser, "--gene_rds", help = "bambu gene RDS.") +parser <- add_argument(parser, "--sample_sheet", help = "Sample sheet CSV.") +parser <- add_argument(parser, "--condition_column", help = "Primary condition column.", default = "condition") +parser <- add_argument(parser, "--covariates", help = "Comma-separated nuisance covariates.") +parser <- add_argument(parser, "--reference_level", help = "Reference level for the condition column.") +parser <- add_argument(parser, "--out_dir", help = "Output directory.", default = "de_analysis") +argv <- parse_args(parser) + +arg_missing <- function(value) { + if (is.null(value) || length(value) == 0 || all(is.na(value))) { + return(TRUE) + } + if (is.character(value)) { + return(all(!nzchar(value))) + } + FALSE +} + +required_args <- c("transcript_rds", "gene_rds", "sample_sheet") +missing_args <- required_args[vapply(required_args, function(arg_name) { + value <- argv[[arg_name]] + arg_missing(value) +}, logical(1))] +if (length(missing_args) > 0) { + stop(sprintf( + "Missing required arguments: %s", + paste(sprintf("--%s", missing_args), collapse = ", ") + )) +} + +dir.create(argv$out_dir, showWarnings = FALSE, recursive = TRUE) + +tx_se <- readRDS(argv$transcript_rds) +gene_se <- readRDS(argv$gene_rds) +sample_df <- read.csv(argv$sample_sheet, check.names = FALSE, stringsAsFactors = FALSE) + +if (!"alias" %in% names(sample_df)) { + stop("Sample sheet must contain an 'alias' column.") +} +if (!(argv$condition_column %in% names(sample_df))) { + stop(sprintf("Sample sheet must contain the '%s' column.", argv$condition_column)) +} + +covariates <- character(0) +if (!arg_missing(argv$covariates)) { + covariates <- trimws(strsplit(argv$covariates, ",", fixed = TRUE)[[1]]) + covariates <- covariates[nzchar(covariates)] +} +missing_covariates <- setdiff(covariates, names(sample_df)) +if (length(missing_covariates) > 0) { + stop(sprintf("Missing covariate columns: %s", paste(missing_covariates, collapse = ", "))) +} + +sample_df <- sample_df[match(colnames(tx_se), sample_df$alias), , drop = FALSE] +if (any(is.na(sample_df$alias))) { + stop("Sample sheet aliases do not match the bambu output sample names.") +} + +condition_values <- unique(sample_df[[argv$condition_column]]) +if (length(condition_values) < 2) { + stop("Differential analysis requires at least two condition levels.") +} + +reference_level <- argv$reference_level +if (arg_missing(reference_level)) { + if ("control" %in% condition_values) { + reference_level <- "control" + } else { + stop("Provide --reference_level when the condition column does not contain 'control'.") + } +} +if (!(reference_level %in% condition_values)) { + stop("The requested reference level is not present in the condition column.") +} + +sample_df[[argv$condition_column]] <- factor(sample_df[[argv$condition_column]]) +for (covariate in covariates) { + sample_df[[covariate]] <- factor(sample_df[[covariate]]) +} + +run_deseq_with_fallback <- function(dds, contrast_name = "unknown") { + tryCatch( + DESeq(dds, quiet = TRUE), + error = function(err) { + if (!grepl( + "all gene-wise dispersion estimates are within 2 orders of magnitude", + conditionMessage(err), + fixed = TRUE + )) { + stop(err) + } + + warning( + "STATISTICAL POWER REDUCED: DESeq2 dispersion estimation failed for ", contrast_name, ".\n", + "This usually indicates:\n", + " 1. Too few replicates (recommend n>=3 per group)\n", + " 2. High biological variability\n", + " 3. Poor data quality\n", + "Falling back to gene-wise dispersion (no information sharing).\n", + "Results will have reduced power and wider confidence intervals." + ) + + dds <- estimateSizeFactors(dds) + dds <- estimateDispersionsGeneEst(dds) + dispersions(dds) <- mcols(dds)$dispGeneEst + + # Write diagnostic file + diag_content <- c( + "DESeq2 Dispersion Estimation Fallback Applied", + "==============================================", + "", + sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")), + sprintf("Contrast: %s", contrast_name), + sprintf("Samples: %d", ncol(dds)), + sprintf("Genes tested: %d", nrow(dds)), + sprintf("Dispersion range: %.3f to %.3f", min(dispersions(dds)), max(dispersions(dds))), + "", + "WHAT HAPPENED:", + " Curve fitting failed. Using gene-wise dispersion estimates.", + "", + "IMPLICATIONS:", + " - No information sharing across genes", + " - Reduced statistical power", + " - Wider confidence intervals", + " - More conservative results (fewer discoveries)", + "", + "LIKELY CAUSES:", + " 1. Too few replicates (recommend n>=3 per group)", + " 2. High biological variability", + " 3. Poor data quality or outlier samples", + "", + "RECOMMENDATIONS:", + " - Add more biological replicates if possible", + " - Check sample quality metrics", + " - Consider filtering low-count genes more stringently" + ) + + diag_file <- file.path(argv$out_dir, sprintf("DESeq2_dispersion_fallback_%s.txt", gsub("[^A-Za-z0-9_-]", "_", contrast_name))) + writeLines(diag_content, diag_file) + + nbinomWaldTest(dds) + } + ) +} + +estimate_dispersions_with_fallback <- function(object, context_label, allow_gene_est = TRUE) { + tryCatch( + estimateDispersions(object), + error = function(err) { + if (!grepl( + "all gene-wise dispersion estimates are within 2 orders of magnitude", + conditionMessage(err), + fixed = TRUE + )) { + stop(err) + } + + message( + context_label, + " dispersion fitting failed; ", + "retrying with fitType='local'." + ) + tryCatch( + estimateDispersions(object, fitType = "local"), + error = function(local_err) { + if (!grepl( + "all gene-wise dispersion estimates are within 2 orders of magnitude", + conditionMessage(local_err), + fixed = TRUE + )) { + stop(local_err) + } + + message( + context_label, + " local-fit dispersion retry failed; ", + "retrying with fitType='mean'." + ) + tryCatch( + estimateDispersions(object, fitType = "mean"), + error = function(mean_err) { + if (!grepl( + "all gene-wise dispersion estimates are within 2 orders of magnitude", + conditionMessage(mean_err), + fixed = TRUE + )) { + stop(mean_err) + } + if (!allow_gene_est) { + stop(mean_err) + } + + message( + context_label, + " mean-fit dispersion retry failed; ", + "falling back to gene-wise dispersion estimates." + ) + object <- estimateDispersionsGeneEst(object) + dispersions(object) <- mcols(object)$dispGeneEst + object + } + ) + } + ) + } + ) +} + +normalise_tsv_value <- function(value) { + if (length(value) == 0 || all(is.na(value))) { + return(NA_character_) + } + if (is.list(value)) { + value <- unlist(value, recursive = TRUE, use.names = FALSE) + } + if (length(value) == 0 || all(is.na(value))) { + return(NA_character_) + } + paste(as.character(value), collapse = ";") +} + +normalise_tsv_df <- function(df) { + as.data.frame( + lapply(df, function(column) { + if (is.list(column)) { + vapply(column, normalise_tsv_value, character(1)) + } else { + column + } + }), + stringsAsFactors = FALSE, + check.names = FALSE + ) +} + +is_recoverable_dexseq_error <- function(message_text) { + grepl( + "all gene-wise dispersion estimates are within 2 orders of magnitude", + message_text, + fixed = TRUE + ) || grepl( + "model matrix is not full rank", + message_text, + fixed = TRUE + ) || grepl( + "replacement has 1 row, data has 0", + message_text, + fixed = TRUE + ) +} + +empty_tsv <- function(columns) { + out <- as.data.frame(matrix(nrow = 0, ncol = length(columns))) + names(out) <- columns + out +} + +write_placeholder_pdf <- function(path, label) { + pdf(path) + plot.new() + text(0.5, 0.5, label, cex = 0.9) + dev.off() +} + +run_deseq2 <- function(count_mat, coldata, target_level, contrast_name) { + design_terms <- c(covariates, argv$condition_column) + design_formula <- as.formula(paste("~", paste(design_terms, collapse = " + "))) + dds <- DESeqDataSetFromMatrix( + countData = round(count_mat), + colData = coldata, + design = design_formula + ) + dds <- run_deseq_with_fallback(dds, contrast_name) + results(dds, contrast = c(argv$condition_column, target_level, reference_level), independentFiltering = TRUE) +} + +run_dexseq <- function(tx_counts, tx_meta, coldata, active_covariates = covariates) { + coldata$sample <- factor(coldata$alias) + coldata[[argv$condition_column]] <- factor(coldata[[argv$condition_column]]) + for (covariate in active_covariates) { + coldata[[covariate]] <- factor(coldata[[covariate]]) + } + + covariate_exon_terms <- if (length(active_covariates) > 0) { + paste0(active_covariates, ":exon") + } else { + character(0) + } + design_terms <- c("sample", "exon", covariate_exon_terms, paste0(argv$condition_column, ":exon")) + reduced_terms <- c("sample", "exon", covariate_exon_terms) + full_formula <- as.formula(paste("~", paste(design_terms, collapse = " + "))) + reduced_formula <- as.formula(paste("~", paste(reduced_terms, collapse = " + "))) + + tryCatch({ + dxd <- DEXSeqDataSet( + countData = round(tx_counts), + sampleData = as.data.frame(coldata), + design = full_formula, + featureID = tx_meta$TXNAME, + groupID = tx_meta$GENEID + ) + dxd <- estimateSizeFactors(dxd) + dxd <- estimate_dispersions_with_fallback(dxd, "DEXSeq", allow_gene_est = TRUE) + dxd <- testForDEU(dxd, reducedModel = reduced_formula) + dxd <- estimateExonFoldChanges(dxd, fitExpToVar = argv$condition_column) + dxr <- DEXSeqResults(dxd, independentFiltering = FALSE) + list(dxd = dxd, dxr = dxr) + }, error = function(err) { + if (length(active_covariates) == 0 || !grepl( + "model matrix is not full rank", + conditionMessage(err), + fixed = TRUE + )) { + stop(err) + } + + dropped_covariate <- tail(active_covariates, 1) + kept_covariates <- head(active_covariates, -1) + message( + "DEXSeq design was not full rank with covariate '", + dropped_covariate, + "'; retrying without it." + ) + run_dexseq(tx_counts, tx_meta, coldata, kept_covariates) + }) +} + +tx_meta <- as.data.frame(rowData(tx_se)) +if (!"TXNAME" %in% names(tx_meta)) { + tx_meta$TXNAME <- rownames(tx_se) +} +if (!"GENEID" %in% names(tx_meta)) { + stop("Transcript rowData must contain GENEID for DEXSeq.") +} + +gene_meta <- as.data.frame(rowData(gene_se)) +if (!"GENEID" %in% names(gene_meta)) { + gene_meta$GENEID <- rownames(gene_se) +} + +targets <- setdiff(as.character(condition_values), reference_level) + +# Initialize QC statistics collector +de_qc_stats <- list() +de_qc_stats$timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S") +de_qc_stats$total_samples <- nrow(sample_df) +de_qc_stats$condition_column <- argv$condition_column +de_qc_stats$reference_level <- reference_level +de_qc_stats$covariates <- if (length(covariates) > 0) covariates else "none" +de_qc_stats$num_contrasts <- length(targets) +de_qc_stats$contrasts <- list() + +# Check sample sizes and warn if underpowered +n_per_group <- table(sample_df[[argv$condition_column]]) +de_qc_stats$samples_per_group <- as.list(n_per_group) + +sample_size_warnings <- c() +if (any(n_per_group < 3)) { + warning( + "WARNING: Some condition groups have fewer than 3 replicates.\n", + "Recommended minimum for DGE: n=3 per group\n", + "Current sample sizes: ", paste(names(n_per_group), "=", n_per_group, collapse=", "), "\n", + "Results may have reduced statistical power." + ) + sample_size_warnings <- c(sample_size_warnings, "Some groups have n<3 (recommended minimum)") +} + +if (any(n_per_group < 2)) { + stop("ERROR: Some condition groups have fewer than 2 replicates. Cannot perform statistical testing.") +} + +de_qc_stats$sample_size_warnings <- if (length(sample_size_warnings) > 0) sample_size_warnings else "none" + +# Multiple testing warning +if (length(targets) > 1) { + fwer <- (1 - (1-0.05)^length(targets)) * 100 + mt_warning <- sprintf( + "Multiple contrasts tested (%d). Per-contrast FDR < 0.05 yields family-wise error rate of ~%.1f%%", + length(targets), fwer + ) + message("WARNING: ", mt_warning) + de_qc_stats$multiple_testing_note <- mt_warning + + mt_content <- c( + "Multiple Testing Across Contrasts", + "==================================", + "", + sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")), + sprintf("Number of contrasts tested: %d", length(targets)), + sprintf("Contrasts: %s", paste(sprintf("%s vs %s", targets, reference_level), collapse=", ")), + "", + "PER-CONTRAST FDR THRESHOLD: 0.05", + sprintf("FAMILY-WISE ERROR RATE: ~%.1f%%", fwer), + "", + "WHAT THIS MEANS:", + " Each contrast uses FDR < 0.05 independently.", + " When testing multiple contrasts, the overall false positive rate increases.", + sprintf(" Expected: %.1f%% chance of at least one false positive across all contrasts", fwer), + "", + "RECOMMENDATIONS:", + " 1. Use stricter per-contrast threshold:", + sprintf(" Bonferroni correction: 0.05 / %d = %.4f", length(targets), 0.05/length(targets)), + " 2. Focus on pre-specified contrasts of interest", + " 3. Treat results as exploratory and validate key findings", + " 4. Consider using hierarchical testing procedures", + "", + "INTERPRETATION:", + " - Results passing FDR < 0.05 in each contrast are discoveries for that contrast", + " - But the overall false discovery burden is higher than 5%", + " - Prioritize genes significant across multiple contrasts", + " - Validate top findings experimentally" + ) + writeLines(mt_content, file.path(argv$out_dir, "MULTIPLE_TESTING_WARNING.txt")) +} + +for (target_level in targets) { + contrast_name <- sprintf("%s_%s_vs_%s", argv$condition_column, target_level, reference_level) + contrast_dir <- file.path(argv$out_dir, contrast_name) + dir.create(contrast_dir, showWarnings = FALSE, recursive = TRUE) + + keep_samples <- sample_df[[argv$condition_column]] %in% c(reference_level, target_level) + contrast_samples <- droplevels(sample_df[keep_samples, , drop = FALSE]) + contrast_samples[[argv$condition_column]] <- relevel( + factor(contrast_samples[[argv$condition_column]]), + ref = reference_level + ) + + # Collect per-contrast QC stats + contrast_qc <- list() + contrast_qc$name <- contrast_name + contrast_qc$target_level <- target_level + contrast_qc$reference_level <- reference_level + contrast_qc$n_samples <- nrow(contrast_samples) + contrast_qc$n_target <- sum(contrast_samples[[argv$condition_column]] == target_level) + contrast_qc$n_reference <- sum(contrast_samples[[argv$condition_column]] == reference_level) + + # DTU power warning + if (nrow(contrast_samples) < 6) { + dtu_warning <- sprintf( + "DTU analysis may be underpowered (n=%d, recommend n>=6 with >=3 per group)", + nrow(contrast_samples) + ) + warning(dtu_warning) + contrast_qc$dtu_power_warning <- dtu_warning + } + + gene_counts <- assays(gene_se)$counts[, contrast_samples$alias, drop = FALSE] + tx_counts <- assays(tx_se)$counts[, contrast_samples$alias, drop = FALSE] + + contrast_qc$genes_tested <- nrow(gene_counts) + contrast_qc$transcripts_tested <- nrow(tx_counts) + + dge_res <- as.data.frame(run_deseq2(gene_counts, contrast_samples, target_level, contrast_name)) + dge_res$GENEID <- rownames(dge_res) + dge_res <- merge(gene_meta, dge_res, by = "GENEID", all.y = TRUE, sort = FALSE) + dge_res <- normalise_tsv_df(dge_res) + + # Collect DGE statistics + contrast_qc$dge_total_genes <- nrow(dge_res) + contrast_qc$dge_significant_fdr05 <- sum(dge_res$padj < 0.05, na.rm = TRUE) + contrast_qc$dge_significant_fdr01 <- sum(dge_res$padj < 0.01, na.rm = TRUE) + contrast_qc$dge_upregulated <- sum(dge_res$padj < 0.05 & dge_res$log2FoldChange > 0, na.rm = TRUE) + contrast_qc$dge_downregulated <- sum(dge_res$padj < 0.05 & dge_res$log2FoldChange < 0, na.rm = TRUE) + + write.table( + dge_res[order(dge_res$padj), ], + file = file.path(contrast_dir, "results_dge.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + + pdf(file.path(contrast_dir, "results_dge.pdf")) + dds_plot <- DESeqDataSetFromMatrix( + countData = round(gene_counts), + colData = contrast_samples, + design = as.formula(paste("~", paste(c(covariates, argv$condition_column), collapse = " + "))) + ) + dds_plot <- run_deseq_with_fallback(dds_plot, contrast_name) + plotMA(results(dds_plot, contrast = c(argv$condition_column, target_level, reference_level), independentFiltering = TRUE)) + dev.off() + + dex_res <- tryCatch( + run_dexseq(tx_counts, tx_meta, contrast_samples), + error = function(err) { + message_text <- conditionMessage(err) + if (!is_recoverable_dexseq_error(message_text)) { + stop(err) + } + + warning( + "DEXSeq failed for contrast ", target_level, " vs ", reference_level, "\n", + "Error: ", message_text + ) + + # Write explicit failure report + failure_content <- c( + "DTU Analysis Failed", + "===================", + "", + sprintf("Timestamp: %s", format(Sys.time(), "%Y-%m-%d %H:%M:%S")), + sprintf("Contrast: %s vs %s", target_level, reference_level), + sprintf("Samples: %d (%d %s, %d %s)", + nrow(contrast_samples), + sum(contrast_samples[[argv$condition_column]] == target_level), target_level, + sum(contrast_samples[[argv$condition_column]] == reference_level), reference_level), + sprintf("Transcripts: %d", nrow(tx_counts)), + "", + "ERROR MESSAGE:", + sprintf(" %s", message_text), + "", + "DTU RESULTS CANNOT BE INTERPRETED", + "", + "This failure is likely due to:", + " 1. Insufficient samples (need >=3 per group, recommend >=6 total for DTU)", + " 2. Too few transcripts with sufficient counts", + " 3. Design matrix not full rank (covariate confounding)", + " 4. Extreme count distributions", + "", + "RECOMMENDATIONS:", + " - Use gene-level DGE results (less power required)", + " - Add more biological replicates", + " - Filter transcripts more stringently", + " - Simplify experimental design (remove problematic covariates)", + "", + "NOTE: Empty DTU result files indicate analysis failure, not 'no DTU detected'" + ) + writeLines(failure_content, file.path(contrast_dir, "DTU_ANALYSIS_FAILED.txt")) + + NULL + } + ) + + if (is.null(dex_res)) { + dex_df <- empty_tsv(c( + "featureID", + "groupID", + "log2fold", + "pvalue", + "padj", + "exonBaseMean" + )) + tx_dtu <- dex_df + gene_dtu <- empty_tsv(c("GENEID", "qval")) + write_placeholder_pdf( + file.path(contrast_dir, "results_dtu.pdf"), + "DEXSeq did not converge for this contrast.\nSee DTU_ANALYSIS_FAILED.txt for details." + ) + contrast_qc$dtu_status <- "FAILED" + contrast_qc$dtu_significant_transcripts <- 0 + contrast_qc$dtu_significant_genes <- 0 + } else { + dxr <- dex_res$dxr + dxd <- dex_res$dxd + dex_df <- as.data.frame(dxr) + dex_df <- normalise_tsv_df(dex_df) + tx_dtu <- dex_df[, intersect( + c("featureID", "groupID", "log2fold", "pvalue", "padj", "exonBaseMean"), + names(dex_df) + ), drop = FALSE] + tx_dtu <- normalise_tsv_df(tx_dtu) + + gene_q <- perGeneQValue(dxr) + gene_dtu <- data.frame( + GENEID = names(gene_q), + qval = unname(gene_q), + row.names = NULL + ) + + # Collect DTU statistics + contrast_qc$dtu_status <- "SUCCESS" + contrast_qc$dtu_significant_transcripts <- sum(tx_dtu$padj < 0.05, na.rm = TRUE) + contrast_qc$dtu_significant_genes <- sum(gene_dtu$qval < 0.05, na.rm = TRUE) + + pdf(file.path(contrast_dir, "results_dtu.pdf")) + plotMA(dxr, cex = 0.8, alpha = 0.05) + plotDispEsts(dxd) + dev.off() + } + + write.table( + dex_df, + file = file.path(contrast_dir, "results_dexseq.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + + write.table( + tx_dtu[order(tx_dtu$padj), ], + file = file.path(contrast_dir, "results_dtu_transcript.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + + write.table( + gene_dtu[order(gene_dtu$qval), ], + file = file.path(contrast_dir, "results_dtu_gene.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + + write.table( + contrast_samples, + file = file.path(contrast_dir, "samples_used.tsv"), + sep = "\t", + quote = FALSE, + row.names = FALSE + ) + + # Write per-contrast QC summary + contrast_qc_summary <- c( + sprintf("Contrast QC Summary: %s", contrast_name), + paste(rep("=", 50), collapse = ""), + "", + "Sample Information:", + sprintf(" Target level (%s): %d samples", target_level, contrast_qc$n_target), + sprintf(" Reference level (%s): %d samples", reference_level, contrast_qc$n_reference), + sprintf(" Total samples: %d", contrast_qc$n_samples), + "", + "DGE Results:", + sprintf(" Genes tested: %d", contrast_qc$genes_tested), + sprintf(" Significant (FDR < 0.05): %d", contrast_qc$dge_significant_fdr05), + sprintf(" Significant (FDR < 0.01): %d", contrast_qc$dge_significant_fdr01), + sprintf(" Upregulated: %d", contrast_qc$dge_upregulated), + sprintf(" Downregulated: %d", contrast_qc$dge_downregulated), + "", + "DTU Results:", + sprintf(" Status: %s", contrast_qc$dtu_status), + sprintf(" Transcripts tested: %d", contrast_qc$transcripts_tested), + if (contrast_qc$dtu_status == "SUCCESS") { + c( + sprintf(" Significant transcripts (FDR < 0.05): %d", contrast_qc$dtu_significant_transcripts), + sprintf(" Genes with DTU (q < 0.05): %d", contrast_qc$dtu_significant_genes) + ) + } else { + " See DTU_ANALYSIS_FAILED.txt for details" + }, + if (!is.null(contrast_qc$dtu_power_warning)) paste0(" WARNING: ", contrast_qc$dtu_power_warning) else NULL, + "" + ) + writeLines(contrast_qc_summary, file.path(contrast_dir, "contrast_qc_summary.txt")) + + # Add to overall QC stats + de_qc_stats$contrasts[[contrast_name]] <- contrast_qc +} + +# Write overall DE/DTU QC statistics as JSON for HTML report +write_json( + de_qc_stats, + file.path(argv$out_dir, "de_qc_stats.json"), + pretty = TRUE, + auto_unbox = TRUE +) + +# Write human-readable overall summary +overall_summary <- c( + "Differential Expression/Usage Analysis Summary", + paste(rep("=", 50), collapse = ""), + "", + sprintf("Timestamp: %s", de_qc_stats$timestamp), + sprintf("Total samples: %d", de_qc_stats$total_samples), + sprintf("Condition column: %s", de_qc_stats$condition_column), + sprintf("Reference level: %s", de_qc_stats$reference_level), + sprintf("Covariates: %s", paste(de_qc_stats$covariates, collapse = ", ")), + "", + "Sample Sizes:", + sapply(names(de_qc_stats$samples_per_group), function(grp) { + sprintf(" %s: %d samples", grp, de_qc_stats$samples_per_group[[grp]]) + }), + if (de_qc_stats$sample_size_warnings != "none") paste0(" WARNING: ", de_qc_stats$sample_size_warnings) else NULL, + "", + sprintf("Number of contrasts tested: %d", de_qc_stats$num_contrasts), + if (!is.null(de_qc_stats$multiple_testing_note)) paste0(" NOTE: ", de_qc_stats$multiple_testing_note) else NULL, + "", + "Per-Contrast Results:", + sapply(names(de_qc_stats$contrasts), function(cname) { + cqc <- de_qc_stats$contrasts[[cname]] + c( + "", + sprintf(" %s:", cname), + sprintf(" Samples: %d (%d vs %d)", cqc$n_samples, cqc$n_target, cqc$n_reference), + sprintf(" DGE significant: %d genes (FDR<0.05)", cqc$dge_significant_fdr05), + sprintf(" DTU status: %s", cqc$dtu_status), + if (cqc$dtu_status == "SUCCESS") sprintf(" DTU significant: %d genes", cqc$dtu_significant_genes) else NULL + ) + }), + "", + "For detailed per-contrast statistics, see:", + " - /contrast_qc_summary.txt", + " - /results_dge.tsv", + " - /results_dtu_gene.tsv", + "" +) +writeLines(unlist(overall_summary), file.path(argv$out_dir, "de_overall_summary.txt")) + +# Save session info for reproducibility +writeLines(capture.output(sessionInfo()), file.path(argv$out_dir, "session_info.txt")) +message("QC statistics written to de_qc_stats.json and de_overall_summary.txt") +message("Session info saved for reproducibility") diff --git a/bin/workflow_glue/check_experiment_design.py b/bin/workflow_glue/check_experiment_design.py new file mode 100644 index 0000000..03b3854 --- /dev/null +++ b/bin/workflow_glue/check_experiment_design.py @@ -0,0 +1,126 @@ +"""Validate DE/DTU sample sheet settings.""" + +from collections import Counter +import csv +import sys + +from .util import get_named_logger, wf_parser # noqa: ABS101 + + +def _split_covariates(value): + if not value: + return [] + return [part.strip() for part in value.split(",") if part.strip()] + + +def main(args): + """Validate sample sheet content for DE/DTU.""" + logger = get_named_logger("checkDesign") + covariates = _split_covariates(args.covariates) + with open(args.sample_sheet, "r", newline="") as handle: + reader = csv.DictReader(handle) + if reader.fieldnames is None: + sys.exit("Sample sheet is empty.") + fieldnames = set(reader.fieldnames) + required = {"alias", args.condition_column} + missing = sorted(required - fieldnames) + if missing: + sys.exit( + "Sample sheet is missing required columns: " + + ", ".join(missing) + ) + + missing_covariates = [name for name in covariates if name not in fieldnames] + if missing_covariates: + sys.exit( + "Sample sheet is missing requested covariate columns: " + + ", ".join(missing_covariates) + ) + + aliases = [] + levels = Counter() + for row in reader: + alias = row.get("alias", "").strip() + if not alias: + sys.exit("Sample sheet contains a row with an empty alias value.") + aliases.append(alias) + + value = row.get(args.condition_column, "").strip() + if not value: + sys.exit( + f"Sample sheet contains a row with an empty " + f"'{args.condition_column}' value." + ) + levels[value] += 1 + + for covariate in covariates: + if not row.get(covariate, "").strip(): + sys.exit( + f"Sample sheet contains an empty value in covariate column " + f"'{covariate}'." + ) + + duplicate_aliases = [ + alias for alias, count in Counter(aliases).items() if count > 1 + ] + if duplicate_aliases: + sys.exit( + "Sample sheet aliases must be unique. Duplicates: " + + ", ".join(sorted(duplicate_aliases)) + ) + + if len(levels) < 2: + sys.exit( + f"The condition column '{args.condition_column}' must contain at " + "least two levels." + ) + + reference_level = args.reference_level + if not reference_level: + if "control" in levels: + reference_level = "control" + else: + sys.exit( + "Provide --reference_level when the condition column does not " + "contain 'control'." + ) + + if reference_level not in levels: + sys.exit( + f"Reference level '{reference_level}' was absent from " + f"'{args.condition_column}'." + ) + + underpowered = [level for level, count in levels.items() if count < 2] + if underpowered: + sys.exit( + "Each condition level must contain at least two samples. " + "Levels with too few samples: " + + ", ".join(sorted(underpowered)) + ) + + logger.info( + "Validated sample sheet %s using condition column '%s' with reference '%s'.", + args.sample_sheet, + args.condition_column, + reference_level, + ) + + +def argparser(): + """Argument parser for the validation entry point.""" + parser = wf_parser("check_experiment_design") + parser.add_argument("--sample_sheet", required=True, help="Sample sheet CSV.") + parser.add_argument( + "--condition_column", default="condition", + help="Primary biological variable column." + ) + parser.add_argument( + "--covariates", default=None, + help="Comma-separated nuisance covariates." + ) + parser.add_argument( + "--reference_level", default=None, + help="Reference level of the primary condition column." + ) + return parser diff --git a/bin/workflow_glue/check_sample_sheet_condition.py b/bin/workflow_glue/check_sample_sheet_condition.py deleted file mode 100644 index afe0ae4..0000000 --- a/bin/workflow_glue/check_sample_sheet_condition.py +++ /dev/null @@ -1,43 +0,0 @@ -#!/usr/bin/env python -"""Check if a sample sheet is valid.""" -from collections import Counter -import csv -import sys - -from .util import get_named_logger, wf_parser # noqa: ABS101 - - -def main(args): - """Run the entry point.""" - logger = get_named_logger("checkSheetCondition") - with open(args.sample_sheet, "r") as f: - csv_reader = csv.DictReader(f) - conditions_count = Counter() - for row in csv_reader: - if "condition" in row: - conditions_count[row['condition']] += 1 - else: - sys.exit( - "Sample sheet has no condition column " - "which is required for the " - "differential expression subworkflow.") - if len(conditions_count.keys()) != 2: - sys.exit( - "There must be only two unique conditions " - "in the condition column of the sample sheet.") - if "control" not in conditions_count: - sys.exit( - "One of the condition types must be control, " - "to indicate which samples to use as the reference.") - if any(v < 2 for v in conditions_count.values()): - sys.exit( - "There must be at least 2 repeats for each " - "condition indicated in the sample sheet.") - logger.info(f"Checked sample sheet for condition column {args.sample_sheet}.") - - -def argparser(): - """Argument parser for entrypoint.""" - parser = wf_parser("check_sample_sheet_condition") - parser.add_argument("sample_sheet", help="Sample sheet to check") - return parser diff --git a/bin/workflow_glue/de_plots.py b/bin/workflow_glue/de_plots.py deleted file mode 100644 index bc41845..0000000 --- a/bin/workflow_glue/de_plots.py +++ /dev/null @@ -1,287 +0,0 @@ -#!/usr/bin/env python -"""Create de report section.""" -import json -import os - -from dominate.tags import h5, p -from dominate.util import raw -from ezcharts import scatterplot -from ezcharts.components.ezchart import EZChart -from ezcharts.layout.snippets import DataTable -from natsort import natsorted -import numpy as np -import pandas as pd - - -def flagstats_df(flagstats_reports): - """Flag stats alignment dataframe.""" - flagstats_dic = {} - for flagstat in flagstats_reports.iterdir(): - with open(flagstat, "r") as f: - data = json.load(f) - data = data["QC-passed reads"] - flagstats = [ - 'mapped', 'primary mapped', 'secondary', 'supplementary'] - per_sample_flagstats = {key: data.get(key) for key in flagstats} - sample = os.path.basename(flagstat).split(".")[0] - flagstats_dic[sample] = per_sample_flagstats - alignment_summary_df = pd.DataFrame(flagstats_dic) - alignment_summary_df = alignment_summary_df[ - natsorted(alignment_summary_df.columns)] - alignment_summary_df.index = [ - "Total Read Mappings", - "Primary", "Secondary", - "Supplementary"] - alignment_summary_df.index.name = "Statistic" - return alignment_summary_df - - -def dexseq_section(dexseq_file, tr_id_to_gene_name, tr_id_to_gene_id, pval_thresh): - """Add gene isoforms table and plot.""" - h5("Differential Isoform usage") - p("""Table showing gene isoforms, ranked by adjusted - p-value, from the DEXSeq analysis. Information shown includes the log2 fold - change between experimental conditions, the log-scaled transcript - abundance and the false discovery corrected p-value (FDR - Benjamini-Hochberg) . - This table has not been filtered - for genes that satisfy statistical or magnitudinal thresholds""") - - dexseq_results = pd.read_csv(dexseq_file, sep='\t') - dexseq_results.index.name = "gene_id:transcript_id" - - # Replace any occurrences of stringtie-generated MSTRG gene ids with - # reference gene_ids. - dexseq_results.index = dexseq_results.index.map( - lambda ge_tr: str( # lookup gene_id from transcript_id [1] - f"{tr_id_to_gene_id.get(ge_tr.split(':')[1])}: {str(ge_tr.split(':')[1])}") - ) - - # Add gene name column. - dexseq_results.insert(0, "gene_name", dexseq_results.index.map( - lambda x: tr_id_to_gene_name.get(x.split(':')[1]))) - - DataTable.from_pandas( - dexseq_results.sort_values(by='pvalue', ascending=True), use_index=True) - - p( - """The figure below presents the MA plot from the DEXSeq analysis. - M is the log2 ratio of isoform transcript abundance between conditions. - A is the log2 transformed mean abundance value. - Transcripts that satisfy the logFC and FDR-corrected - (False discovery rate - Benjamini-Hochberg) p-value - thresholds defined are shaded as 'Up-' or 'Down-' regulated.""") - - dexseq_results['direction'] = 'not_sig' - - dexseq_results.loc[ - (dexseq_results["Log2FC"] > 0) & (dexseq_results['pvalue'] < pval_thresh), - 'direction'] = 'up' - - dexseq_results.loc[ - (dexseq_results["Log2FC"] <= 0) & (dexseq_results['pvalue'] < pval_thresh), - 'direction'] = 'down' - - plot = scatterplot( - data=dexseq_results, x='Log2MeanExon', y='Log2FC', hue='direction', - palette=['#E32636', '#7E8896', '#0A22DE'], - hue_order=['up', 'down', 'not_sig'], marker='circle') - plot._fig.xaxis.axis_label = "A (log2 transformed mean exon read counts)" - plot._fig.yaxis.axis_label = "M (log2 transformed differential abundance)" - plot.legend = dict(orient='horizontal', top=30) - plot._fig.title = "Average copy per million (CPM) vs Log-fold change (LFC)" - EZChart(plot) - - -def dtu_section(dtu_file, txid_to_gene_name): - """Plot dtu section.""" - dtu_results = pd.read_csv(dtu_file, sep='\t') - dtu_results["gene_name"] = dtu_results["txID"].apply( - lambda x: txid_to_gene_name.get(x)) - - dtu_pvals = dtu_results.sort_values(by='gene', ascending=True) - raw("""Table showing gene and transcript identifiers - and their FDR-corrected (False discovery rate - Benjamini-Hochberg) probabilities - for the genes and their isoforms that have been - identified as showing DTU using the R packages DEXSeq and StageR. - This list has been shortened requiring that both gene and transcript - must satisfy the p-value - threshold""") - DataTable.from_pandas(dtu_results.loc[dtu_pvals.index], use_index=False) - - raw("""View dtu_plots.pdf file to see plots of differential isoform usage""") - - -def dge_section(df, pval_thresh): - """Create DGE table and MA plot.""" - h5("Differential gene expression") - df[['logFC', 'logCPM', 'F']] = df[ - ['logFC', 'logCPM', 'F']].round(2) - - p("""Table showing the genes from the edgeR analysis. - Information shown includes the log2 fold change between - experimental conditions, the log-scaled counts per million measure of abundance - and the FDR-corrected p-value (False discovery rate - Benjamini-Hochberg). - This table has not been - filtered for genes that satisfy statistical or magnitudinal thresholds""") - - df = df.sort_values('FDR', ascending=True) - df.index.name = 'gene_id' - DataTable.from_pandas(df, use_index=True) - - h5("Results of the edgeR Analysis.") - - p("""This plot visualises differences in measurements between the - two experimental conditions. M is the log2 ratio of gene expression - calculated between the conditions. - A is a log2 transformed mean expression value. - The figure below presents the MA figure from this edgeR analysis. - Genes that satisfy the logFC and FDR-corrected - (False discovery rate - Benjamini-Hochberg) p-value thresholds - defined are shaded as 'Up-' or 'Down-' regulated. - """) - df['sig'] = None - df.loc[(df["logFC"] > 0) & (df['PValue'] < pval_thresh), 'sig'] = 'up' - df.loc[(df["logFC"] <= 0) & (df['PValue'] < pval_thresh), 'sig'] = 'down' - df.loc[(df["PValue"] >= pval_thresh), 'sig'] = 'not_sig' - - plot = scatterplot( - data=df, x='logCPM', y='logFC', hue='sig', - palette=['#E32636', '#7E8896', '#0A22DE'], - hue_order=['up', 'not_sig', 'down'], marker='circle') - plot._fig.x_range.start = 10 - plot._fig.xaxis.axis_label = "Average log CPM" - plot._fig.yaxis.axis_label = "Log-fold change" - plot.legend = dict(orient='horizontal', top=30) - # Should opacity of the symbols be lowered? - plot._fig.title = "Average copy per million (CPM) vs Log-fold change (LFC)" - EZChart(plot) - - -def salmon_table(salmon_counts): - """Create salmon counts summary table.""" - salmon_counts = pd.read_csv(salmon_counts, sep='\t') - salmon_counts.set_index("Reference", drop=True, append=False, inplace=True) - salmon_size_top = salmon_counts.sum(axis=1).sort_values(ascending=False) - salmon_counts = salmon_counts.applymap(np.int64) - h5("Transcripts Per Million") - p("""Table showing the annotated Transcripts Per Million - identified by Minimap2 mapping and Salmon transcript - detection. Displaying the top 100 transcripts with the highest - number of mapped reads""") - - salmon_counts = salmon_counts[sorted(salmon_counts.columns)] - DataTable.from_pandas( - salmon_counts.loc[salmon_size_top.index].head(n=100), use_index=True) - - -def get_translations(gtf): - """Create gene_and transcript id mappings. - - Annotation can be stringtie-generated (GTF) or from the input - reference annotation (GTF or GFF3) and the various attributes can differ - """ - with open(gtf) as fh: - txid_to_gene_name = {} - gid_to_gene_name = {} - tx_id_to_gene_id = {} - - def get_feature(row, feature): - return row.split(feature)[1].split( - ";")[0].replace('=', '').replace("\"", "").strip() - - for gff_entry in fh: - # Process transcripts features only - if gff_entry.startswith("#") or gff_entry.split('\t')[2] != 'transcript': - continue - # Different gtf/gff formats contain different attributes - # and different formating (eg. gene_name="xyz" or gene_name "xyz") - gene_name = gene_id = transcript_id = 'unknown' - - if 'ref_gene_id' in gff_entry: - # Favour ref_gene_id over gene_id. The latter can be multi-locus merged - # genes from stringtie - gene_id = get_feature(gff_entry, 'ref_gene_id') - elif 'gene_id' in gff_entry: - gene_id = get_feature(gff_entry, 'gene_id') - else: - gene_id = get_feature(gff_entry, 'gene') - - if 'transcript_id' in gff_entry: - transcript_id = get_feature(gff_entry, 'transcript_id') - - if 'gene_name' in gff_entry: - gene_name = get_feature(gff_entry, 'gene_name') - else: - # Fallback to gene_id if gene_name is not present - gene_name = gene_id - - txid_to_gene_name[transcript_id] = gene_name - tx_id_to_gene_id[transcript_id] = gene_id - gid_to_gene_name[gene_id] = gene_name - return txid_to_gene_name, tx_id_to_gene_id, gid_to_gene_name - - -def de_section( - annotation, dge, dexseq, dtu, - tpm, report, filtered, unfiltered, - gene_counts, flagstats_dir, pval_threshold=0.01): - """Differential expression sections.""" - with (report.add_section("Differential expression", "DE")): - - p("""This section shows differential gene expression - and differential isoform usage. Salmon was used to - assign reads to individual annotated isoforms defined by - the GTF-format annotation. - These counts were used to perform a statistical analysis to identify - the genes and isoforms that show differences in abundance between - the experimental conditions. - Any novel genes or transcripts that do not have relevant gene or transcript IDs - are prefixed with MSTRG for use in differential expression analysis. - Find the full sequences of any transcripts in the - final_non_redundant_transcriptome.fasta file. - """) - alignment_summary_df = flagstats_df(flagstats_dir) - h5("Alignment summary stats") - DataTable.from_pandas(alignment_summary_df, use_index=True) - salmon_table(tpm) - - # Get translations for adding gene names to tables - ( - txid_to_gene_name, txid_to_gene_id, gid_to_gene_name - ) = get_translations(annotation) - - # Add gene names columns to counts files and write out - # for publishing to user dir. - df_dge = pd.read_csv(dge, sep='\t') - df_dge.insert(0, 'gene_name', df_dge.index.map( - lambda x: gid_to_gene_name.get(x))) - df_dge.to_csv('results_dge.tsv', index=True, index_label="gene_id", sep="\t") - - df_gene_counts = pd.read_csv(gene_counts, sep='\t') - df_gene_counts.insert( - 0, 'gene_name', df_gene_counts.index.map( - lambda x: gid_to_gene_name.get(x))) - df_gene_counts.to_csv( - 'all_gene_counts.tsv', index=True, index_label="gene_id", sep="\t") - - df_filtered = pd.read_csv(filtered, sep='\t') - df_filtered.insert(1, "gene_name", df_filtered.gene_id.map( - lambda x: gid_to_gene_name.get(x))) - df_filtered.to_csv( - 'filtered_transcript_counts_with_genes.tsv', index=False, sep='\t') - - df_unfiltered = pd.read_csv(unfiltered, sep='\t') - df_unfiltered.insert(1, "gene_name", df_unfiltered.gene_id.map( - lambda x: gid_to_gene_name.get(x))) - df_unfiltered.to_csv( - 'unfiltered_transcript_counts_with_genes.tsv', index=False, sep='\t') - - df_tpm = pd.read_csv(tpm, sep='\t') - df_tpm.insert(1, "gene_name", df_tpm.Reference.map( - lambda x: txid_to_gene_name.get(x))) - df_tpm.to_csv("unfiltered_tpm_transcript_counts.tsv", index=False, sep='\t') - - # Add tables to report - dge_section(df_dge, pval_threshold) - dexseq_section(dexseq, txid_to_gene_name, txid_to_gene_id, pval_threshold) - dtu_section(dtu, txid_to_gene_name) diff --git a/bin/workflow_glue/generate_pychopper_stats.py b/bin/workflow_glue/generate_pychopper_stats.py deleted file mode 100644 index 8955986..0000000 --- a/bin/workflow_glue/generate_pychopper_stats.py +++ /dev/null @@ -1,39 +0,0 @@ -#!/usr/bin/env python -"""Generate CSV of pychopper stats.""" - -# -*- coding: utf-8 -*- - -import os - -import pandas as pd - -from .util import wf_parser # noqa: ABS101 - - -def argparser(): - """Argument parser for entrypoint.""" - parser = wf_parser("generate_pychopper_stats") - parser.add_argument("--data", required=True, help="") - parser.add_argument("--output_dir", required=True, help="") - - return parser - - -def generate_pychopper_stats(tsv, output): - """Make CSV of pychopper stats.""" - classified_path = os.path.join(output, "pychopper_stats.csv") - df = pd.read_csv(tsv, sep="\t", index_col="Name") - classified = df.loc[df["Category"] == "Classification"]\ - .copy().reset_index().rename(columns={'Name': 'Classification'}) - classified["Percentage"] = \ - 100 * classified["Value"] / classified["Value"].sum() - tuning = df.loc[df["Category"] == "AutotuneSample"]\ - .copy().reset_index().rename(columns={'Name': 'Filter'}) - tuning.to_csv(classified_path) - - -def main(args): - """Run entry point.""" - assert os.path.isfile(args.data) - assert os.path.isdir(args.output_dir) - generate_pychopper_stats(tsv=args.data, output=args.output_dir) diff --git a/bin/workflow_glue/merge_count_tsvs.py b/bin/workflow_glue/merge_count_tsvs.py deleted file mode 100644 index 59bd6fa..0000000 --- a/bin/workflow_glue/merge_count_tsvs.py +++ /dev/null @@ -1,64 +0,0 @@ -#!/usr/bin/env python -"""Merge salmon output count files.""" - -from functools import reduce - -import numpy as np -import pandas as pd - -from .util import wf_parser # noqa: ABS101 - - -def argparser(): - """Argument parser for entrypoint.""" - parser = wf_parser("merge_count_tsvs") - parser.add_argument( - '-j', metavar='join', help="Join type (outer).", default="outer") - parser.add_argument( - '-f', metavar='field', - help="Join on this field (Reference).", default="Reference") - parser.add_argument( - '-o', metavar='out_tsv', - help="Output tsv (merge_tsvs.tsv).", default="merge_tsvs.tsv") - parser.add_argument( - '-z', action="store_true", - help="Fill NA values with zero.", default=False) - parser.add_argument( - '-tpm', type=bool, default=False, - help="TPM instead of counts") - parser.add_argument( - '-tsvs', metavar='input_tsvs', nargs='*', - help="Input tab separated files.") - - return parser - - -def main(args): - """Run entry point.""" - dfs = {x: pd.read_csv(x, sep="\t") for x in args.tsvs} - - ndfs = [] - for x, df in dfs.items(): - # Transform counts to integers: - if args.tpm: - df = df.rename(columns={'TPM': 'Count', 'Name': 'Reference'}) - else: - df = df.rename(columns={'NumReads': 'Count', 'Name': 'Reference'}) - df.Count = np.array(df.Count, dtype=int) - # Take only non-zero counts: - df = df[df.Count > 0] - df = df[["Reference", "Count"]] - df = df.sort_values(by=["Count"], ascending=False) - name = x.split('.')[0] - df = df.rename(columns={'Count': name}) - ndfs.append(df) - dfs = ndfs - - df_merged = reduce(lambda left, right: pd.merge( - left, right, on=args.f, how=args.j), dfs) - if args.z: - df_merged = df_merged.fillna(0) - df_merged = df_merged.sort_index(axis=1) - df_merged = df_merged.sort_index(axis=0) - - df_merged.to_csv(args.o, sep="\t", index=False) diff --git a/bin/workflow_glue/models/custom.py b/bin/workflow_glue/models/custom.py new file mode 100644 index 0000000..5fc2235 --- /dev/null +++ b/bin/workflow_glue/models/custom.py @@ -0,0 +1,16 @@ +"""Extended models for the workflow.""" +from dataclasses import dataclass + +from workflow_glue.models.common import CheckResult + + +@dataclass +class CheckResult(CheckResult): + """ + A result of some check the workflow has performed. + + This can be at a sample or workflow level. + """ + + categories = dict( + example_check="Example check category") diff --git a/bin/workflow_glue/parse_gffcompare.py b/bin/workflow_glue/parse_gffcompare.py deleted file mode 100644 index a71630a..0000000 --- a/bin/workflow_glue/parse_gffcompare.py +++ /dev/null @@ -1,264 +0,0 @@ -#!/usr/bin/env python -"""Make report tables and data for plotting.""" -import os -from pathlib import Path - -import numpy as np -import pandas as pd -from .util import get_named_logger, wf_parser # noqa: ABS101 - - -def argparser(): - """Argument parser for entrypoint.""" - parser = wf_parser("Parse gffcompare") - parser.add_argument( - '--sample_id', help="Sample ID", required=True) - parser.add_argument( - '--gffcompare_dir', - help="The gffcompare output directory", - required=False, - type=Path) - parser.add_argument( - '--isoform_table_out', - help="Output path for per-isoform table", - type=Path) - parser.add_argument( - '--tracking', - help="gffcompare tracking file", - type=Path) - parser.add_argument( - "--annotation", - required=False, - default=None, help="Reference annotation GFF file") - - return parser - - -def _parse_stat_line(sl): - """Parse a stats line.""" - res = {} - tmp = sl.split(':')[1].split('|') - res['sensitivity'] = float(tmp[0].strip()) - res['precision'] = float(tmp[1].strip()) - return res - - -def _parse_matching_line(line): - """Parse a matching line.""" - tmp = line.split(':')[1].strip() - return int(tmp) - - -def _parse_mn_line(line): - """Parse a miss or novel line.""" - res = {} - tmp = line.split(':')[1].strip() - tmp = tmp.split('/') - res['value'] = int(tmp[0]) - tmp = tmp[1].split('(') - res['value_total'] = int(tmp[0].strip()) - res['percent'] = float(tmp[1].split('%)')[0]) - return res - - -def _parse_total_line(line): - """Parse a total line.""" - res = {} - tmp = line.split(':')[1].strip() - tmp = tmp.split('in') - res['transcripts'] = int(tmp[0].strip()) - tmp = tmp[1].split('loci') - res['loci'] = int(tmp[0].strip()) - tmp = int(tmp[1].split('(')[1].split(' ')[0]) - res['me_transcripts'] = tmp - return res - - -def parse_gffcmp_stats(gffcompare_stats, sample_id, outpath): - """Parse a gffcompare stats file. - - Gffcompare stats file - :param gffcompare_stats: Path to the gffcompare stats file. - :returns: Return as tuple of dataframes containing: - perfromance statistics, match statistics, miss statistics, - novel statistics, total statistics. - :rtype: tuple - """ - performance = [] - missed = [] - novel = [] - total = [] - - with open(gffcompare_stats, 'r') as fh: - for line in fh: - line = line.strip() - if len(line) == 0: - continue - - # Parse totals: - if line.startswith('# Query mRNAs'): - r = _parse_total_line(line) - total.append([r['loci'], 'loci', 'query']) - total.append([r['transcripts'], 'transcripts', 'query']) - total.append([r['me_transcripts'], 'multexonic', 'query']) - if line.startswith('# Reference mRNAs '): - r = _parse_total_line(line) - total.append([r['loci'], 'loci', 'reference']) - total.append([r['transcripts'], 'transcripts', 'reference']) - total.append([r['me_transcripts'], 'multexonic', 'reference']) - - # Parse basic statistics: - if line.startswith('Base level'): - st = _parse_stat_line(line) - performance.append((st['sensitivity'], 'Sensitivity', 'Base')) - performance.append((st['precision'], 'Precision', 'Base')) - if line.startswith('Exon level'): - st = _parse_stat_line(line) - performance.append((st['sensitivity'], 'Sensitivity', 'Exon')) - performance.append((st['precision'], 'Precision', 'Exon')) - if line.startswith('Intron level'): - st = _parse_stat_line(line) - performance.append((st['sensitivity'], 'Sensitivity', 'Intron')) - performance.append((st['precision'], 'Precision', 'Intron')) - if line.startswith('Intron chain level'): - st = _parse_stat_line(line) - performance.append((st['sensitivity'], 'Sensitivity', 'Intron_chain')) - performance.append((st['precision'], 'Precision', 'Intron_chain')) - if line.startswith('Transcript level'): - st = _parse_stat_line(line) - performance.append((st['sensitivity'], 'Sensitivity', 'Transcript')) - performance.append((st['precision'], 'Precision', 'Transcript')) - if line.startswith('Locus level'): - st = _parse_stat_line(line) - performance.append((st['sensitivity'], 'Sensitivity', 'Locus')) - performance.append((st['precision'], 'Precision', 'Locus')) - - # Parse missing statistics: - if line.startswith('Missed exons'): - r = _parse_mn_line(line) - missed.append((r['value'], 'Missed', 'Exons')) - missed.append((r['value_total'], 'total', 'Exons')) - missed.append((r['percent'], 'Percent', 'Exons')) - if line.startswith('Missed introns'): - r = _parse_mn_line(line) - missed.append((r['value'], 'Missed', 'Introns')) - missed.append((r['value_total'], 'total', 'Introns')) - missed.append((r['percent'], 'Percent', 'Introns')) - if line.startswith('Missed loci'): - r = _parse_mn_line(line) - missed.append((r['value'], 'Missed', 'Loci')) - missed.append((r['value_total'], 'total', 'Loci')) - missed.append((r['percent'], 'Percent', 'Loci')) - - # Parse novel statistics: - if line.startswith('Novel exons'): - r = _parse_mn_line(line) - novel.append((r['value'], 'Novel', 'Exons')) - novel.append((r['value_total'], 'Total', 'Exons')) - novel.append((r['percent'], 'Percent_novel', 'Exons')) - if line.startswith('Novel introns'): - r = _parse_mn_line(line) - novel.append((r['value'], 'Novel', 'Introns')) - novel.append((r['value_total'], 'Total', 'Introns')) - novel.append((r['percent'], 'Percent_novel', 'Introns')) - if line.startswith('Novel loci'): - r = _parse_mn_line(line) - novel.append((r['value'], 'Novel', 'Loci')) - novel.append((r['value_total'], 'Total', 'Loci')) - novel.append((r['percent'], 'Percent_novel', 'Loci')) - - def write_records(records, fn): - pd.DataFrame.from_records(records, columns=['counts', 'type', 'source']) \ - .to_csv(outpath / fn, sep='\t') - - write_records(total, 'Totals.tsv') - write_records(missed, 'Missed.tsv') - write_records(performance, 'Performance.tsv') - write_records(novel, 'Novel.tsv') - - -def tracking_summary(tracking_file, output_dir, annotations=None): - """Write per transcript class gffcompare tracking files.""" - tracking_headings = [ - "query_transfrag_id", "query_locus_id", "ref_gene_id", - "class", "details"] - nice_names = { - '=': 'complete', 'c': 'contained', 'k': 'containment', - 'm': 'retained', 'n': 'retained (partial)', 'j': 'multi', - 'e': 'single', 'o': 'overlap', 's': 'opposite', - 'x': 'exonic', 'i': 'intron', 'y': 'contains', 'p': 'runon', - 'r': 'repeat', 'u': 'unknown'} - - if os.path.exists(annotations): - tracking = pd.read_csv( - tracking_file, sep="\t", names=tracking_headings[1:], - index_col=0) - - df = ( - pd.DataFrame(tracking['class'].value_counts()) - .reset_index() - # Reminder when updating pandas > 2 - # .rename(columns={'count': 'Count'}) - .rename(columns={'index': 'class', 'class': 'Count'}) - ) - - df['Percent'] = round(df['Count'] * 100 / df['Count'].sum(), 2) - df['description'] = [nice_names[x] for x in df['class']] - - df = df.sort_values('Count', ascending=True) - df.to_csv(output_dir / 'tracking_summary.tsv', sep='\t') - - else: - logger = get_named_logger('trackingSum') - logger.info("Skipping classification summary as no annotation provided.") - - -def make_isoform_table(gffcompare_dir, sample_id, outpath): - """Make an isoform summary table.""" - try: - tmap_file = next(gffcompare_dir.glob('*.tmap')) - except StopIteration: - raise ValueError("Cannot find .tmap file in {}".format(gffcompare_dir)) - dtypes = { - 'ref_gene_id': str, - 'ref_id': str, - 'class_code': str, - 'qry_id': str, - 'num_exons': np.uint16, - 'cov': np.uint32, - 'len': np.uint32 - } - df = pd.read_csv( - tmap_file, sep='\t+', - index_col=None, - usecols=list(dtypes.keys()), - dtype=dtypes) - - if df.empty: # No transcripts. Write a header only result file - df = pd.DataFrame( - columns=list(dtypes.keys()) + ['sample_id', 'parent gene iso num']) - df.to_csv(f'{sample_id}_transcripts_table.tsv', sep='\t', index=False) - else: - df = df.assign(sample_id=sample_id) - - # Make a column of number of isoforms in parent gene - gb = df.groupby(['ref_gene_id']).count() - gb.rename(columns={'ref_id': 'num_isoforms'}, inplace=True) - - df['parent gene iso num'] = df.apply( - lambda x: gb.loc[(x.ref_gene_id), 'num_isoforms'], axis=1) - - # Unclassified transcripts should not be lumped together - df.loc[df.class_code == 'u', 'parent gene iso num'] = None - - df.to_csv(outpath, sep='\t', index=False) - - -def main(args): - """Entry point.""" - if args.gffcompare_dir: # TODO: should this every be optional? - stats = args.gffcompare_dir / 'str_merged.stats' - parse_gffcmp_stats(stats, args.sample_id, args.gffcompare_dir) - make_isoform_table(args.gffcompare_dir, args.sample_id, args.isoform_table_out) - tracking_summary( - args.tracking, args.gffcompare_dir, args.annotation) diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py index 14692f4..e77a126 100644 --- a/bin/workflow_glue/report.py +++ b/bin/workflow_glue/report.py @@ -1,403 +1,682 @@ -#!/usr/bin/env python -"""Create workflow report.""" +"""Create workflow report for wf-transcriptomes.""" import json from pathlib import Path -import pickle -from bokeh.models import HoverTool, Range1d -from bokeh.models.tickers import AdaptiveTicker -from dominate.tags import li, ul +from dominate.tags import div, h3, p, pre, strong from dominate.util import raw -import ezcharts as ezc from ezcharts.components import fastcat -from ezcharts.components.ezchart import EZChart -from ezcharts.components.reports.labs import LabsReport -from ezcharts.components.theme import LAB_head_resources -from ezcharts.layout.snippets import DataTable, Grid, Tabs -from ezcharts.plots.categorical import barplot -from ezcharts.util import get_named_logger +from ezcharts.components.reports import labs +from ezcharts.layout.snippets import Tabs +from ezcharts.layout.snippets.table import DataTable import pandas as pd -from . import de_plots # noqa: ABS101 -from .util import wf_parser # noqa: ABS101 + +from .util import get_named_logger, wf_parser # noqa: ABS101 -def argparser(): - """Argument parser for entrypoint.""" - parser = wf_parser("report") - parser.add_argument("--report", help="Report output file") - parser.add_argument( - "--metadata", default='metadata.json', required=True, - help="sample metadata") - parser.add_argument( - "--stats", nargs='+', - help="Fastcat per-read stats, ordered as per entries in --metadata.") - parser.add_argument( - "--versions", required=True, - help="directory containing CSVs containing name,version.") - parser.add_argument( - "--params", default=None, required=True, - help="A JSON file containing the workflow parameter key/values") - parser.add_argument( - "--wf_version", default='unknown', - help="version of the executed workflow") - parser.add_argument( - "--alignment_stats", required=False, default=None, type=Path, - help="TSV summary file of alignment statistics") - parser.add_argument( - "--gff_annotation", required=False, type=Path, - help="transcriptome annotation gff file") - parser.add_argument( - "--gffcompare_dir", required=False, default=None, type=Path, - help="gffcompare outout dir") - parser.add_argument( - "--pychop_report", required=False, default=None, type=Path, - help="TSV summary file of pychopper statistics") - parser.add_argument( - "--isoform_table", required=False, type=Path, - help="Path to directory of TSV files with isoform summaries") - parser.add_argument( - "--isoform_table_nrows", required=False, type=int, default=5000, - help="Maximum rows to display in isoforms table") - parser.add_argument( - "--transcriptome_summary", required=False, type=Path, - help="Path to dir containing transcriptome summary results ") - parser.add_argument( - "--de_report", required=False, type=Path, default=None, - help="Differential expression report optional") - parser.add_argument( - "--de_stats", required=False, type=Path, default=None, - help="Differential expression report optional") - parser.add_argument( - "--pval_threshold", required=False, type=float, default=0.01, - help=( - "pvalue theshold for inclusion of differentially expressed genes" - " transcripts in plots ")) - - return parser +def _find_table(directory, pattern): + matches = sorted(Path(directory).glob(pattern)) + return matches[0] if matches else None -def gff_compare_plots(report, gffcompare_outdirs): - """Create various sections and plots in a WfReport. - - :param gffcompare_outdirs: List of output directories from run_gffcompare - :return: None - - """ - # Plot overview panel: - with report.add_section("Annotation summary", "Annotation"): - - raw("""The following plots summarize some of the output from - gffcompare - """) - - raw("""