Add configurable attribute_type parameter for featureCounts

.github/workflows/test-PR.yml

@@ -8,22 +8,22 @@ 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
           cd tests/test_counts/reads/
           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
 

.github/workflows/test_push.yml

@@ -11,22 +11,22 @@ 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
           cd tests/test_counts/reads/
           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()

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

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"
 

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

envs/gffread.yaml

@@ -0,0 +1,6 @@
+name: gffread
+channels:
+  - conda-forge
+  - bioconda
+dependencies:
+  - gffread=0.12.7

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";

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"""

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} \\

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} \\

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()