diff --git a/.github/workflows/test-PR.yml b/.github/workflows/test-PR.yml index 339cd4d..cfc9050 100644 --- a/.github/workflows/test-PR.yml +++ b/.github/workflows/test-PR.yml @@ -8,12 +8,12 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v3 - + - name: Setup Python for FTP server uses: actions/setup-python@v4 with: python-version: '3.x' - + - name: Install pyftpdlib and start FTP server run: | pip install pyftpdlib @@ -21,9 +21,9 @@ jobs: python ftp_server.py & echo $! > $GITHUB_WORKSPACE/ftp_server.pid sleep 3 - curl ftp://localhost:21/ || echo "FTP server connection test failed" + curl ftp://localhost:2121/ || echo "FTP server connection test failed" echo "FTP server started with PID: $(cat $GITHUB_WORKSPACE/ftp_server.pid)" - + - name: Install Miniforge and Snakemake shell: bash -l {0} run: | @@ -40,27 +40,31 @@ jobs: run: | source "${HOME}/conda/etc/profile.d/conda.sh" conda activate snakemake - + echo "Running dry-run test..." - snakemake --cores 10 --use-conda --dry-run --config \ + snakemake --cores 4 --use-conda --dry-run --config \ processing_dir="tests/test_counts/test_out/processing" \ results_dir="tests/test_counts/test_out/results" \ benchmark_dir="tests/test_counts/test_out/benchmarks" \ reference_fasta="tests/test_counts/genome/random_genome.fa" \ gtf_file="tests/test_counts/genome/annotation.gtf" \ samples_csv="test_samples_long.csv" \ - cleanup_processing='True' + cleanup_processing='True' \ + cores_align=2 cores_coverage=2 cores_fastp=2 cores_featurecounts=2 \ + cores_flagstat=2 cores_nanopore_align=2 max_cores=2 echo "Running full test with FTP..." - snakemake --cores 10 --use-conda --config \ + snakemake --cores 4 --use-conda --config \ processing_dir="tests/test_counts/test_out/processing" \ results_dir="tests/test_counts/test_out/results" \ benchmark_dir="tests/test_counts/test_out/benchmarks" \ reference_fasta="tests/test_counts/genome/random_genome.fa" \ gtf_file="tests/test_counts/genome/annotation.gtf" \ samples_csv="test_samples_long.csv" \ - cleanup_processing='True' - + cleanup_processing='True' \ + cores_align=2 cores_coverage=2 cores_fastp=2 cores_featurecounts=2 \ + cores_flagstat=2 cores_nanopore_align=2 max_cores=2 + ls tests/test_counts/test_out/results -l head tests/test_counts/test_out/results/SAMPLE_02/SAMPLE_02_counts_single_all.txt diff --git a/.github/workflows/test_push.yml b/.github/workflows/test_push.yml index 784efc6..6b1fcd9 100644 --- a/.github/workflows/test_push.yml +++ b/.github/workflows/test_push.yml @@ -11,12 +11,12 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v3 - + - name: Setup Python for FTP server uses: actions/setup-python@v4 with: python-version: '3.x' - + - name: Install pyftpdlib and start FTP server run: | pip install pyftpdlib @@ -24,9 +24,9 @@ jobs: python ftp_server.py & echo $! > $GITHUB_WORKSPACE/ftp_server.pid sleep 3 - curl ftp://localhost:21/ || echo "FTP server connection test failed" + curl ftp://localhost:2121/ || echo "FTP server connection test failed" echo "FTP server started with PID: $(cat $GITHUB_WORKSPACE/ftp_server.pid)" - + - name: Install Miniforge and Snakemake shell: bash -l {0} run: | @@ -43,17 +43,19 @@ jobs: run: | source "${HOME}/conda/etc/profile.d/conda.sh" conda activate snakemake - + echo "Running dry-run test..." - snakemake --cores 10 --use-conda --dry-run --config \ + snakemake --cores 4 --use-conda --dry-run --config \ processing_dir="tests/test_counts/test_out/processing" \ results_dir="tests/test_counts/test_out/results" \ benchmark_dir="tests/test_counts/test_out/benchmarks" \ reference_fasta="tests/test_counts/genome/random_genome.fa" \ gtf_file="tests/test_counts/genome/annotation.gtf" \ - samples_csv="test_samples_long.csv" \ - cleanup_processing='True' - + samples_csv="test_samples_ftp_and_local.csv" \ + cleanup_processing='True' \ + cores_align=2 cores_coverage=2 cores_fastp=2 cores_featurecounts=2 \ + cores_flagstat=2 cores_nanopore_align=2 max_cores=2 + - name: Stop FTP server if: always() diff --git a/README.md b/README.md index c70897b..6d74385 100644 --- a/README.md +++ b/README.md @@ -122,6 +122,18 @@ Bowtie2 index automatically, so you only need to provide the reference FASTA (with a `.fa` extension) and annotation file. The pipeline will exit with an error if the FASTA does not use the required `.fa` suffix. +**Annotation Format Support:** +The pipeline accepts both GTF and GFF annotation files. When a GFF file is provided, +it is implicitly converted to GTF format during processing. +```bash +gffread in.gff3 -T -o out.gtf +``` +Users should ensure that +the attributes and feature types specified in the configuration (e.g., `feature_type`, +`attribute_type`) are present in the original file and will be preserved in the +GFF-to-GTF conversion. This is particularly important for custom annotations or +non-standard attribute names. + ### 4. Run the pipeline ```bash @@ -130,6 +142,7 @@ snakemake -n # Run the pipeline snakemake --cores 8 +``` ### Running the pipeline with the provided Singularity container @@ -156,7 +169,6 @@ snakemake --use-singularity \ Additional details about the container, including the software stack, are available in the companion repository: [mtinti/myRna-seq-docker](https://github.com/mtinti/myRna-seq-docker). -``` ## Testing the Pipeline @@ -169,14 +181,13 @@ results_dir="tests/test_counts/new_branch/results" \ benchmark_dir="tests/test_counts/new_branch/benchmarks" \ reference_fasta="tests/test_counts/genome/random_genome.fa" \ gtf_file="tests/test_counts/genome/annotation.gtf" \ -samples_csv="test_samples_local.csv" -``` +samples_csv="test_samples_local.csv"``` Make sure the FASTA you reference ends with `.fa`; other extensions will cause the workflow to abort during validation. Expected output from the test dataset: -``` + Features --------------------------- gene1: 10 reads total diff --git a/Snakefile b/Snakefile index 9ffdd8a..952c34c 100644 --- a/Snakefile +++ b/Snakefile @@ -106,8 +106,22 @@ except Exception as e: if "processing_genome_index" not in config: fasta_name = os.path.splitext(os.path.basename(config["reference_fasta"]))[0] config["processing_genome_index"] = os.path.join(config['processing_dir'], 'reference', fasta_name) + if "processing_annotation_source" not in config: + config["processing_annotation_source"] = os.path.join( + config['processing_dir'], 'reference', os.path.basename(config["gtf_file"]) + ) if "processing_gtf_file" not in config: - config["processing_gtf_file"] = os.path.join(config['processing_dir'], 'reference', os.path.basename(config["gtf_file"])) + gtf_basename = os.path.basename(config["gtf_file"]) + annotation_stem, annotation_ext = os.path.splitext(gtf_basename) + if annotation_stem.endswith(".gff"): + annotation_stem = os.path.splitext(annotation_stem)[0] + if annotation_ext.lower() in (".gff", ".gff3"): + derived_gtf = f"{annotation_stem}.gtf" + else: + derived_gtf = gtf_basename + config["processing_gtf_file"] = os.path.join( + config['processing_dir'], 'reference', derived_gtf + ) if "processing_reference_fasta" not in config: config["processing_reference_fasta"] = os.path.join( config['processing_dir'], 'reference', os.path.basename(config["reference_fasta"]) or "reference.fa" @@ -203,6 +217,9 @@ onstart: # Include common functions include: "rules/common.smk" +# Prepare annotation for downstream steps (convert GFF->GTF, validate feature/attribute) +include: "rules/processing/prepare_annotation.smk" + # Build Bowtie2 indexes from the staged reference FASTA include: "rules/processing/build_bowtie_index.smk" diff --git a/config.yaml b/config.yaml index 009e03f..0b1e09e 100644 --- a/config.yaml +++ b/config.yaml @@ -23,6 +23,7 @@ cores_fastp: 8 cores_default: 1 max_cores: 8 feature_type: "CDS" +attribute_type: "gene_id" # Attribute used for grouping features (e.g., gene_id, transcript_id, ID) # Nanopore-specific settings nanopore_preset: "map-ont" # Minimap2 preset for nanopore diff --git a/envs/gffread.yaml b/envs/gffread.yaml new file mode 100644 index 0000000..df2a187 --- /dev/null +++ b/envs/gffread.yaml @@ -0,0 +1,6 @@ +name: gffread +channels: + - conda-forge + - bioconda +dependencies: + - gffread=0.12.7 diff --git a/my.gtf b/my.gtf new file mode 100644 index 0000000..4a8bd26 --- /dev/null +++ b/my.gtf @@ -0,0 +1,9 @@ +test_genome Prediction transcript 11 500 . - . transcript_id "gene1:mRNA"; gene_id "gene1" +test_genome Prediction exon 11 500 . - . transcript_id "gene1:mRNA"; gene_id "gene1"; +test_genome Prediction CDS 11 500 . - 0 transcript_id "gene1:mRNA"; gene_id "gene1"; +test_genome Prediction transcript 2001 2500 . - . transcript_id "gene2:mRNA"; gene_id "gene1" +test_genome Prediction exon 2001 2500 . - . transcript_id "gene2:mRNA"; gene_id "gene1"; +test_genome Prediction CDS 2001 2500 . - 0 transcript_id "gene2:mRNA"; gene_id "gene1"; +test_genome Prediction transcript 4001 4500 . - . transcript_id "gene3:mRNA"; gene_id "gene1" +test_genome Prediction exon 4001 4500 . - . transcript_id "gene3:mRNA"; gene_id "gene1"; +test_genome Prediction CDS 4001 4500 . - 0 transcript_id "gene3:mRNA"; gene_id "gene1"; diff --git a/rules/__pycache__/sample_utils.cpython-313.pyc b/rules/__pycache__/sample_utils.cpython-313.pyc new file mode 100644 index 0000000..533d8ca Binary files /dev/null and b/rules/__pycache__/sample_utils.cpython-313.pyc differ diff --git a/rules/common.smk b/rules/common.smk index 8ffe26d..c8c4aa0 100644 --- a/rules/common.smk +++ b/rules/common.smk @@ -37,6 +37,11 @@ def get_bowtie2_index_files(): """Return the expected Bowtie2 index files within the processing directory.""" return [f"{config['processing_genome_index']}{ext}" for ext in BOWTIE2_INDEX_SUFFIXES] + +def get_annotation_gtf(wildcards=None): + """Return the staged annotation path used across the pipeline.""" + return config["processing_gtf_file"] + # Helper function to get results directory path def get_results_path(*args): """Get a path within the results directory specified in config""" diff --git a/rules/counting/paired_counts.smk b/rules/counting/paired_counts.smk index 1d4ac0f..b337f68 100644 --- a/rules/counting/paired_counts.smk +++ b/rules/counting/paired_counts.smk @@ -39,6 +39,7 @@ rule featurecounts_paired_all: params: gtf = config["processing_gtf_file"], feature_type = config["feature_type"], + attribute_type = config["attribute_type"], extra = "-p -B -C -M -O --countReadPairs" # Paired-end specific options # -B requireBothEndsMapped # -p isPairedEnd @@ -81,7 +82,7 @@ rule featurecounts_paired_all: featureCounts {params.extra} \\ -T {threads} \\ -t {params.feature_type} \\ - -g gene_id \\ + -g {params.attribute_type} \\ -a {params.gtf} \\ -o {output.counts} \\ {input.bam} \\ @@ -121,7 +122,8 @@ rule featurecounts_paired_unique: flag = get_processing_path("{run_tag}/featurecounts_paired_complete.flag") params: gtf = config["processing_gtf_file"], - feature_type = config["feature_type"], + feature_type = config["feature_type"], + attribute_type = config["attribute_type"], extra = "-p -B -C -O --countReadPairs", # Paired-end specific options min_mapping_quality = config.get("min_mapping_quality", 2) log: @@ -161,7 +163,7 @@ rule featurecounts_paired_unique: featureCounts {params.extra} \\ -T {threads} \\ -t {params.feature_type} \\ - -g gene_id \\ + -g {params.attribute_type} \\ -Q {params.min_mapping_quality} \\ -a {params.gtf} \\ -o {output.counts} \\ diff --git a/rules/counting/single_counts.smk b/rules/counting/single_counts.smk index fe9b040..2aeefe7 100644 --- a/rules/counting/single_counts.smk +++ b/rules/counting/single_counts.smk @@ -38,7 +38,8 @@ rule featurecounts_single_all: summary = get_processing_path("{run_tag}/qc/feature_counts/{run_tag}_counts_single_all.txt.summary") params: gtf = config["processing_gtf_file"], - feature_type = config["feature_type"], + feature_type = config["feature_type"], + attribute_type = config["attribute_type"], extra = "-C -O -M" # Single-end specific options (no -p, -B) log: get_processing_path("{run_tag}/logs/featurecounts_single_all.log") @@ -76,7 +77,7 @@ rule featurecounts_single_all: featureCounts {params.extra} \\ -T {threads} \\ -t {params.feature_type} \\ - -g gene_id \\ + -g {params.attribute_type} \\ -a {params.gtf} \\ -o {output.counts} \\ {input.bam} \\ @@ -117,6 +118,7 @@ rule featurecounts_single_unique: params: gtf = config["processing_gtf_file"], feature_type = config["feature_type"], + attribute_type = config["attribute_type"], extra = "-C -O", # Single-end specific options (no -p, -B) min_mapping_quality = config.get("min_mapping_quality", 2) log: @@ -156,7 +158,7 @@ rule featurecounts_single_unique: featureCounts {params.extra} \\ -T {threads} \\ -t {params.feature_type} \\ - -g gene_id \\ + -g {params.attribute_type} \\ -Q {params.min_mapping_quality} \\ -a {params.gtf} \\ -o {output.counts} \\ diff --git a/rules/processing/prepare_annotation.smk b/rules/processing/prepare_annotation.smk new file mode 100644 index 0000000..4c1c1f3 --- /dev/null +++ b/rules/processing/prepare_annotation.smk @@ -0,0 +1,61 @@ +""" +Prepare the annotation for downstream rules. + +* If the user supplies a GFF/GFF3 annotation, convert it to GTF via gffread. +* Validate that the configured feature and attribute exist in the final GTF. +""" + +import os +from pathlib import Path + +from sample_utils import validate_feature_attribute + +ANNOTATION_SOURCE = config["processing_annotation_source"] +ANNOTATION_GTF = config["processing_gtf_file"] +VALIDATION_FLAG = os.path.join(config["processing_dir"], "reference", ".annotation_validated") +ANNOTATION_NEEDS_CONVERSION = ANNOTATION_SOURCE != ANNOTATION_GTF + + +def _validate(gtf_path, feature_type=None, attribute_type=None): + validate_feature_attribute( + gtf_path, + feature_type or config.get("feature_type"), + attribute_type or config.get("attribute_type"), + ) + + +if ANNOTATION_NEEDS_CONVERSION: + rule prepare_annotation: + input: + source=ANNOTATION_SOURCE, + output: + gtf=ANNOTATION_GTF, + flag=VALIDATION_FLAG, + params: + feature_type=lambda wildcards: config.get("feature_type"), + attribute_type=lambda wildcards: config.get("attribute_type"), + threads: + 1 + conda: + "../../envs/gffread.yaml" + singularity: + config.get("singularity_image", "") + run: + shell("mkdir -p $(dirname {output.gtf})") + shell("gffread {input.source} -T -o {output.gtf}") + _validate(output.gtf, params.feature_type, params.attribute_type) + Path(output.flag).touch() +else: + rule prepare_annotation: + input: + gtf=ANNOTATION_GTF, + output: + flag=VALIDATION_FLAG, + params: + feature_type=lambda wildcards: config.get("feature_type"), + attribute_type=lambda wildcards: config.get("attribute_type"), + threads: + 1 + run: + _validate(input.gtf, params.feature_type, params.attribute_type) + Path(output.flag).touch() diff --git a/rules/qc/paired_rnaseq.smk b/rules/qc/paired_rnaseq.smk index a9b8ecb..0d02f4b 100644 --- a/rules/qc/paired_rnaseq.smk +++ b/rules/qc/paired_rnaseq.smk @@ -9,13 +9,13 @@ rule qualimap_rnaseq_paired: input: bam = get_picard_bam, bai = lambda wildcards: f"{get_picard_bam(wildcards)}.bai", - bamqc_flag = get_processing_path("{run_tag}/bamqc_complete.flag") + bamqc_flag = get_processing_path("{run_tag}/bamqc_complete.flag"), + gtf = get_annotation_gtf output: dir = directory(get_processing_path("{run_tag}/qc/{run_tag}_qualimap_rnaseq/")), flag = get_processing_path("{run_tag}/qc_paired_complete.flag") params: - memory = config.get("qualimap_memory", "10G"), - gtf = config["processing_gtf_file"] + memory = config.get("qualimap_memory", "10G") log: stderr = get_processing_path("{run_tag}/logs/qualimap_rnaseq.log"), stdout = get_processing_path("{run_tag}/logs/qualimap_rnaseq.stdout.log") @@ -36,20 +36,20 @@ rule qualimap_rnaseq_paired: # Log start echo "Running qualimap rnaseq (paired-end) on {wildcards.run_tag}" > {log.stderr} echo "Input BAM: {input.bam}" >> {log.stderr} - echo "GTF file: {params.gtf}" >> {log.stderr} + echo "GTF file: {input.gtf}" >> {log.stderr} echo "Output directory: {output.dir}" >> {log.stderr} echo "Java memory: {params.memory}" >> {log.stderr} - + # Check if GTF file exists - if [[ ! -f {params.gtf} ]]; then - echo "ERROR: GTF file {params.gtf} not found" >> {log.stderr} + if [[ ! -f {input.gtf} ]]; then + echo "ERROR: GTF file {input.gtf} not found" >> {log.stderr} exit 1 fi - + # Run qualimap rnaseq with paired-end option qualimap rnaseq --java-mem-size={params.memory} \\ -bam {input.bam} \\ - -gtf {params.gtf} \\ + -gtf {input.gtf} \\ -outdir {output.dir} \\ -outformat HTML \\ -pe \\ diff --git a/rules/qc/single_rnaseq.smk b/rules/qc/single_rnaseq.smk index 2e31a68..fa2891c 100644 --- a/rules/qc/single_rnaseq.smk +++ b/rules/qc/single_rnaseq.smk @@ -9,13 +9,13 @@ rule qualimap_rnaseq_single: input: bam = get_picard_bam, bai = lambda wildcards: f"{get_picard_bam(wildcards)}.bai", - bamqc_flag = get_processing_path("{run_tag}/bamqc_complete.flag") + bamqc_flag = get_processing_path("{run_tag}/bamqc_complete.flag"), + gtf = get_annotation_gtf output: dir = directory(get_processing_path("{run_tag}/qc/{run_tag}_qualimap_rnaseq/")), flag = get_processing_path("{run_tag}/qc_single_complete.flag") params: - memory = config.get("qualimap_memory", "10G"), - gtf = config["processing_gtf_file"] + memory = config.get("qualimap_memory", "10G") log: stderr = get_processing_path("{run_tag}/logs/qualimap_rnaseq.log"), stdout = get_processing_path("{run_tag}/logs/qualimap_rnaseq.stdout.log") @@ -36,20 +36,20 @@ rule qualimap_rnaseq_single: # Log start echo "Running qualimap rnaseq (single-end) on {wildcards.run_tag}" > {log.stderr} echo "Input BAM: {input.bam}" >> {log.stderr} - echo "GTF file: {params.gtf}" >> {log.stderr} + echo "GTF file: {input.gtf}" >> {log.stderr} echo "Output directory: {output.dir}" >> {log.stderr} echo "Java memory: {params.memory}" >> {log.stderr} - + # Check if GTF file exists - if [[ ! -f {params.gtf} ]]; then - echo "ERROR: GTF file {params.gtf} not found" >> {log.stderr} + if [[ ! -f {input.gtf} ]]; then + echo "ERROR: GTF file {input.gtf} not found" >> {log.stderr} exit 1 fi - + # Run qualimap rnaseq without paired-end option qualimap rnaseq --java-mem-size={params.memory} \\ -bam {input.bam} \\ - -gtf {params.gtf} \\ + -gtf {input.gtf} \\ -outdir {output.dir} \\ -outformat HTML \\ > {log.stdout} 2>> {log.stderr} diff --git a/rules/sample_utils.py b/rules/sample_utils.py index f1939ba..8a7d1b8 100644 --- a/rules/sample_utils.py +++ b/rules/sample_utils.py @@ -1,6 +1,7 @@ """ Utility functions for sample handling in the RNA-seq pipeline. """ +import gzip import os import pandas as pd import shutil @@ -35,25 +36,52 @@ def copy_reference_files(config): fasta_stem = os.path.splitext(fasta_stem)[0] target_fasta = os.path.join(target_dir, fasta_basename) target_genome_base = os.path.join(target_dir, fasta_stem) - target_gtf_file = os.path.join(target_dir, os.path.basename(src_gtf_file)) + + annotation_basename = os.path.basename(src_gtf_file) + annotation_stem, annotation_ext = os.path.splitext(annotation_basename) + if annotation_stem.endswith(".gff"): + annotation_stem = os.path.splitext(annotation_stem)[0] + + target_annotation_source = os.path.join(target_dir, annotation_basename) + if annotation_ext.lower() in (".gff", ".gff3"): + target_gtf_file = os.path.join(target_dir, f"{annotation_stem}.gtf") + else: + target_gtf_file = target_annotation_source config["processing_reference_fasta"] = target_fasta config["processing_genome_index"] = target_genome_base + config["processing_annotation_source"] = target_annotation_source config["processing_gtf_file"] = target_gtf_file + # Copy FASTA if it doesn't exist or if source is newer + needs_fasta_copy = False if not os.path.exists(target_fasta): + needs_fasta_copy = True + elif os.path.getmtime(src_fasta) > os.path.getmtime(target_fasta): + needs_fasta_copy = True + print(f"Source FASTA is newer than cached copy, updating...") + + if needs_fasta_copy: print(f"Copying FASTA file from {src_fasta} to {target_fasta}") try: shutil.copy2(src_fasta, target_fasta) except Exception as e: print(f"Warning: Failed to copy FASTA file {src_fasta}: {e}") - if not os.path.exists(target_gtf_file): - print(f"Copying GTF file from {src_gtf_file} to {target_gtf_file}") + # Copy annotation if it doesn't exist or if source is newer + needs_annotation_copy = False + if not os.path.exists(target_annotation_source): + needs_annotation_copy = True + elif os.path.getmtime(src_gtf_file) > os.path.getmtime(target_annotation_source): + needs_annotation_copy = True + print(f"Source annotation is newer than cached copy, updating...") + + if needs_annotation_copy: + print(f"Copying annotation file from {src_gtf_file} to {target_annotation_source}") try: - shutil.copy2(src_gtf_file, target_gtf_file) + shutil.copy2(src_gtf_file, target_annotation_source) except Exception as e: - print(f"Warning: Failed to copy GTF file: {e}") + print(f"Warning: Failed to copy annotation file: {e}") def create_sample_directories(config, sample_name): """Create the necessary directory structure for a specific sample""" @@ -283,6 +311,58 @@ def validate_reference_inputs(config): print(f"Found reference FASTA at {fasta_path} and GTF at {gtf_path}") + +def validate_feature_attribute(gtf_path, feature_type, attribute_key): + """Ensure the requested feature and attribute exist in a GTF/GFF3 file.""" + + if not feature_type or not attribute_key: + raise ValueError( + "ERROR: feature_type and attribute_type must be set in the config to validate the annotation." + ) + + found_feature = False + found_attribute = False + open_func = gzip.open if gtf_path.endswith(".gz") else open + + with open_func(gtf_path, "rt") as handle: + for line in handle: + if not line or line.startswith("#"): + continue + + parts = line.rstrip("\n").split("\t") + if len(parts) < 9 or parts[2] != feature_type: + continue + + found_feature = True + attributes = parts[8] + + for field in attributes.split(";"): + field = field.strip() + if not field: + continue + + if " " in field: + key, _ = field.split(" ", 1) + elif "=" in field: + key, _ = field.split("=", 1) + else: + continue + + if key == attribute_key: + found_attribute = True + break + + if found_attribute: + break + + if not found_feature: + raise ValueError(f"ERROR: Feature type '{feature_type}' not found in {gtf_path}.") + + if not found_attribute: + raise ValueError( + f"ERROR: Attribute '{attribute_key}' not found for feature type '{feature_type}' in {gtf_path}." + ) + def is_sample_completed(config, sample, samples_df): """Check if a sample has already been processed""" try: diff --git a/test_config.yaml b/test_config.yaml index a206ac6..c15bcbc 100644 --- a/test_config.yaml +++ b/test_config.yaml @@ -23,6 +23,7 @@ cores_fastp: 8 cores_default: 1 max_cores: 8 feature_type: "CDS" +attribute_type: "gene_id" # Attribute used for grouping features (e.g., gene_id, transcript_id, ID) # Mark Duplicates parameters remove_duplicates: False # Set to true to remove duplicates, false to only mark them diff --git a/tests/.DS_Store b/tests/.DS_Store new file mode 100644 index 0000000..384b3be Binary files /dev/null and b/tests/.DS_Store differ diff --git a/tests/test_counts/genome/.DS_Store b/tests/test_counts/genome/.DS_Store new file mode 100644 index 0000000..5008ddf Binary files /dev/null and b/tests/test_counts/genome/.DS_Store differ diff --git a/tests/test_counts/genome/test.gff b/tests/test_counts/genome/test.gff new file mode 100644 index 0000000..f9c401d --- /dev/null +++ b/tests/test_counts/genome/test.gff @@ -0,0 +1,15 @@ +##gff-version 3 +##species http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=185431 +##sequence-region test_genome 1 5000 +test_genome Prediction gene 11 500 . - . ID=gene1;description=test_desc gene1;ebi_biotype=protein_coding +test_genome Prediction mRNA 11 500 . - . ID=gene1:mRNA;Parent=gene1;description=test_desc gene1;gene_ebi_biotype=protein_coding +test_genome Prediction exon 11 500 . - . ID=exon_gene1:mRNA-E1;Parent=gene1:mRNA;gene_id=gene1 +test_genome Prediction CDS 11 500 . - 0 ID=gene1:mRNA-p1-CDS1;Parent=gene1:mRNA;gene_id=gene1;protein_source_id=gene1:mRNA-p1 +test_genome Prediction gene 2001 2500 . - . ID=gene2;description=test_desc gene2;ebi_biotype=protein_coding +test_genome Prediction mRNA 2001 2500 . - . ID=gene2:mRNA;Parent=gene1;description=test_desc gene2;gene_ebi_biotype=protein_coding +test_genome Prediction exon 2001 2500 . - . ID=exon_gene2:mRNA-E1;Parent=gene2:mRNA;gene_id=gene2 +test_genome Prediction CDS 2001 2500 . - 0 ID=gene2:mRNA-p1-CDS1;Parent=gene2:mRNA;gene_id=gene2;protein_source_id=gene2:mRNA-p1 +test_genome Prediction gene 4001 4500 . - . ID=gene3;description=test_desc gene3;ebi_biotype=protein_coding +test_genome Prediction mRNA 4001 4500 . - . ID=gene3:mRNA;Parent=gene1;description=test_desc gene3;gene_ebi_biotype=protein_coding +test_genome Prediction exon 4001 4500 . - . ID=exon_gene3:mRNA-E1;Parent=gene3:mRNA;gene_id=gene3 +test_genome Prediction CDS 4001 4500 . - 0 ID=gene3:mRNA-p1-CDS1;Parent=gene3:mRNA;gene_id=gene3;protein_source_id=gene3:mRNA-p1 diff --git a/tests/test_counts/gff_config.yaml b/tests/test_counts/gff_config.yaml new file mode 100644 index 0000000..3a12649 --- /dev/null +++ b/tests/test_counts/gff_config.yaml @@ -0,0 +1,37 @@ +# Test configuration using GFF annotation to exercise conversion +processing_dir: "tests/test_counts/test_out_gff/processing" +results_dir: "tests/test_counts/test_out_gff/results" +benchmark_dir: "tests/test_counts/test_out_gff/benchmarks" + +samples_csv: "test_samples_long.csv" +selected_samples: [] + +reference_fasta: "tests/test_counts/genome/random_genome.fa" +gtf_file: "tests/test_counts/genome/test.gff" + +cores_align: 4 +cores_coverage: 4 +cores_flagstat: 4 +cores_featurecounts: 4 +cores_fastp: 4 +cores_default: 1 +max_cores: 4 +feature_type: "CDS" +attribute_type: "gene_id" + +remove_duplicates: False +qualimap_memory: "2G" + +coverage_bin_size: 50 +coverage_normalize: "RPKM" +min_mapping_quality: 2 + +copy_fastq: False +copy_bam: False +copy_benchmarks: False +copy_logs: False +cleanup_processing: False + +max_resources: + network: 2 + io: 1000