Add RPKM calculation step; harden RPKM script against zero denominators

.github/workflows/test-PR.yml

@@ -67,6 +67,7 @@ jobs:
 
           ls tests/test_counts/test_out/results -l
           head tests/test_counts/test_out/results/SAMPLE_02/SAMPLE_02_counts_single_all.txt
+          head tests/test_counts/test_out/results/SAMPLE_02/SAMPLE_02_rpkm_single_all.txt
 
       - name: Stop FTP server
         if: always()

README.md

@@ -38,6 +38,8 @@ Coverage Tracks (bamCoverage)
     ↓
 Feature Counting (featureCounts)
     ↓
+RPKM Calculation
+    ↓
 Results & Cleanup
 
 * Step skipped for nanopore data

Snakefile

@@ -270,6 +270,7 @@ conditional_modules = [
     ("rules/coverage/single_coverage.smk", HAS_SINGLE_OR_NANO_EFFECTIVE),
     ("rules/counting/paired_counts.smk", HAS_PAIRED_EFFECTIVE),
     ("rules/counting/single_counts.smk", HAS_SINGLE_OR_NANO_EFFECTIVE),
+    ("rules/counting/rpkm.smk", HAS_PAIRED_EFFECTIVE or HAS_SINGLE_OR_NANO_EFFECTIVE),
 ]
 
 for module, should_include in conditional_modules:
@@ -303,8 +304,9 @@ workflow_steps = [
     "07 - BAM Quality Control",
     "08 - Coverage Tracks",
     "09 - Feature Counting",
-    "10 - Benchmark Analysis",
-    "11 - Copy Results"
+    "10 - RPKM Calculation",
+    "11 - Benchmark Analysis",
+    "12 - Copy Results"
 ]
 
 onsuccess:
@@ -322,4 +324,4 @@ onsuccess:
 onerror:
     print("RNA-seq pipeline encountered an error.")
     print("Check the log files in the processing directory for details.")
-    print("You can resume the pipeline with --keep-going to continue from where it left off.")
+    print("You can resume the pipeline with --keep-going to continue from where it left off.")

envs/rpkm.yaml

@@ -0,0 +1,8 @@
+name: rpkm
+channels:
+  - conda-forge
+  - defaults
+dependencies:
+  - python=3.11
+  - pandas
+  - numpy

rules/common.smk

@@ -260,6 +260,13 @@ def get_featurecounts_flag(run_tag):
     else:
         return get_processing_path(f"{run_tag}/featurecounts_single_complete.flag")
 
+def get_rpkm_flag(run_tag):
+    """Get the appropriate RPKM flag based on run tag type"""
+    if is_paired_end(run_tag):
+        return get_processing_path(f"{run_tag}/rpkm_paired_complete.flag")
+    else:
+        return get_processing_path(f"{run_tag}/rpkm_single_complete.flag")
+
 def get_markduplicates_flag(wildcards):
     """Get the mark duplicates flag for a run tag or sample."""
     if isinstance(wildcards, dict):
@@ -274,4 +281,4 @@ def get_markduplicates_flag(wildcards):
 
 def get_sample_cleanup_flag(run_tag):
     """Get the per-sample cleanup flag for a run tag"""
-    return get_results_path(f"{run_tag}/sample_cleanup_complete.txt")
+    return get_results_path(f"{run_tag}/sample_cleanup_complete.txt")

rules/counting/rpkm.smk

@@ -0,0 +1,127 @@
+"""
+Rules for calculating RPKM values from featureCounts outputs.
+Generates RPKM files for both all reads and uniquely mapped reads.
+"""
+import re
+
+
+def is_valid_paired_run_tag(run_tag):
+    if run_tag in SAMPLES_DF.index and run_tag not in RUN_TAGS:
+        return is_paired_end(run_tag)
+    if run_tag in RUN_TAG_SAMPLES:
+        samples = RUN_TAG_SAMPLES[run_tag]
+        if samples:
+            return is_paired_end(samples[0])
+    return False
+
+
+def is_valid_single_run_tag(run_tag):
+    if run_tag in SAMPLES_DF.index and run_tag not in RUN_TAGS:
+        return is_single_end(run_tag) or is_nanopore(run_tag)
+    if run_tag in RUN_TAG_SAMPLES:
+        samples = RUN_TAG_SAMPLES[run_tag]
+        if samples:
+            return is_single_end(samples[0]) or is_nanopore(samples[0])
+    return False
+
+
+PAIRED_RUN_TAGS = [rt for rt in EFFECTIVE_SAMPLES_TO_PROCESS if is_valid_paired_run_tag(rt)]
+SINGLE_RUN_TAGS = [rt for rt in EFFECTIVE_SAMPLES_TO_PROCESS if is_valid_single_run_tag(rt)]
+
+
+rule rpkm_paired_all:
+    wildcard_constraints:
+        run_tag = "|".join([re.escape(rt) for rt in PAIRED_RUN_TAGS]) if PAIRED_RUN_TAGS else "^$"
+    input:
+        counts = get_processing_path("{run_tag}/{run_tag}_counts_paired_all.txt")
+    output:
+        rpkm = get_processing_path("{run_tag}/{run_tag}_rpkm_paired_all.txt")
+    log:
+        get_processing_path("{run_tag}/logs/rpkm_paired_all.log")
+    conda:
+        "../../envs/rpkm.yaml"
+    script:
+        "../scripts/calc_rpkm.py"
+
+
+rule rpkm_paired_unique:
+    wildcard_constraints:
+        run_tag = "|".join([re.escape(rt) for rt in PAIRED_RUN_TAGS]) if PAIRED_RUN_TAGS else "^$"
+    input:
+        counts = get_processing_path("{run_tag}/{run_tag}_counts_paired_unique.txt")
+    output:
+        rpkm = get_processing_path("{run_tag}/{run_tag}_rpkm_paired_unique.txt")
+    log:
+        get_processing_path("{run_tag}/logs/rpkm_paired_unique.log")
+    conda:
+        "../../envs/rpkm.yaml"
+    script:
+        "../scripts/calc_rpkm.py"
+
+
+rule rpkm_paired_complete:
+    wildcard_constraints:
+        run_tag = "|".join([re.escape(rt) for rt in PAIRED_RUN_TAGS]) if PAIRED_RUN_TAGS else "^$"
+    input:
+        rpkm_all = get_processing_path("{run_tag}/{run_tag}_rpkm_paired_all.txt"),
+        rpkm_unique = get_processing_path("{run_tag}/{run_tag}_rpkm_paired_unique.txt")
+    output:
+        flag = get_processing_path("{run_tag}/rpkm_paired_complete.flag")
+    log:
+        get_processing_path("{run_tag}/logs/rpkm_paired_complete.log")
+    shell:
+        """
+        echo "RPKM calculation completed for paired-end sample {wildcards.run_tag}" > {output.flag}
+        echo "Timestamp: $(date)" >> {output.flag}
+        echo "All reads RPKM: {input.rpkm_all}" >> {output.flag}
+        echo "Unique reads RPKM: {input.rpkm_unique}" >> {output.flag}
+        """
+
+
+rule rpkm_single_all:
+    wildcard_constraints:
+        run_tag = "|".join([re.escape(rt) for rt in SINGLE_RUN_TAGS]) if SINGLE_RUN_TAGS else "^$"
+    input:
+        counts = get_processing_path("{run_tag}/{run_tag}_counts_single_all.txt")
+    output:
+        rpkm = get_processing_path("{run_tag}/{run_tag}_rpkm_single_all.txt")
+    log:
+        get_processing_path("{run_tag}/logs/rpkm_single_all.log")
+    conda:
+        "../../envs/rpkm.yaml"
+    script:
+        "../scripts/calc_rpkm.py"
+
+
+rule rpkm_single_unique:
+    wildcard_constraints:
+        run_tag = "|".join([re.escape(rt) for rt in SINGLE_RUN_TAGS]) if SINGLE_RUN_TAGS else "^$"
+    input:
+        counts = get_processing_path("{run_tag}/{run_tag}_counts_single_unique.txt")
+    output:
+        rpkm = get_processing_path("{run_tag}/{run_tag}_rpkm_single_unique.txt")
+    log:
+        get_processing_path("{run_tag}/logs/rpkm_single_unique.log")
+    conda:
+        "../../envs/rpkm.yaml"
+    script:
+        "../scripts/calc_rpkm.py"
+
+
+rule rpkm_single_complete:
+    wildcard_constraints:
+        run_tag = "|".join([re.escape(rt) for rt in SINGLE_RUN_TAGS]) if SINGLE_RUN_TAGS else "^$"
+    input:
+        rpkm_all = get_processing_path("{run_tag}/{run_tag}_rpkm_single_all.txt"),
+        rpkm_unique = get_processing_path("{run_tag}/{run_tag}_rpkm_single_unique.txt")
+    output:
+        flag = get_processing_path("{run_tag}/rpkm_single_complete.flag")
+    log:
+        get_processing_path("{run_tag}/logs/rpkm_single_complete.log")
+    shell:
+        """
+        echo "RPKM calculation completed for single-end sample {wildcards.run_tag}" > {output.flag}
+        echo "Timestamp: $(date)" >> {output.flag}
+        echo "All reads RPKM: {input.rpkm_all}" >> {output.flag}
+        echo "Unique reads RPKM: {input.rpkm_unique}" >> {output.flag}
+        """

rules/results.smk

@@ -14,6 +14,7 @@ rule copy_results:
         # Required inputs from previous steps - use helper functions to get type-specific flags
         qc_complete = lambda wildcards: get_qc_flag(wildcards.run_tag),
         featurecounts_complete = lambda wildcards: get_featurecounts_flag(wildcards.run_tag),
+        rpkm_complete = lambda wildcards: get_rpkm_flag(wildcards.run_tag),
         coverage_complete = lambda wildcards: get_coverage_flag(wildcards.run_tag),
         benchmark_summary = get_processing_path("{run_tag}/benchmarks_summary.txt")
     output:
@@ -93,6 +94,10 @@ rule copy_results:
         # Copy count files
         echo "Copying count files" >> {log}
         cp {params.proc_dir}/{wildcards.run_tag}_counts_*.txt {params.results_sample_dir}/ 2>> {log}
+
+        # Copy RPKM files
+        echo "Copying RPKM files" >> {log}
+        cp {params.proc_dir}/{wildcards.run_tag}_rpkm_*.txt {params.results_sample_dir}/ 2>> {log}
 
         # Copy benchmark summary file
         echo "Copying benchmark summary" >> {log}
@@ -181,6 +186,7 @@ rule copy_results:
         echo "- QC results" >> {output.complete}
         echo "- Coverage tracks (bigWig)" >> {output.complete}
         echo "- Feature counts" >> {output.complete}
+        echo "- RPKM values" >> {output.complete}
         echo "- Benchmark summary" >> {output.complete}
 
         if [ "{params.copy_bam}" = "true" ]; then
@@ -401,4 +407,4 @@ rule copy_results_all:
                 echo "Reference directory does not exist: {params.reference_dir}" >> {log}
             fi
         fi
-        """
+        """

rules/sample_utils.py

@@ -381,6 +381,10 @@ def is_sample_completed(config, sample, samples_df):
             # Check if counts files exist
             paired_counts_all = os.path.join(config['results_dir'], effective_name, f"{effective_name}_counts_paired_all.txt")
             paired_counts_unique = os.path.join(config['results_dir'], effective_name, f"{effective_name}_counts_paired_unique.txt")
+
+            # Check if RPKM files exist
+            paired_rpkm_all = os.path.join(config['results_dir'], effective_name, f"{effective_name}_rpkm_paired_all.txt")
+            paired_rpkm_unique = os.path.join(config['results_dir'], effective_name, f"{effective_name}_rpkm_paired_unique.txt")
 
             # Check if coverage files exist
             all_bw = os.path.join(config['results_dir'], effective_name, f"{effective_name}_all.bw")
@@ -390,13 +394,19 @@ def is_sample_completed(config, sample, samples_df):
             final_results_exist = (
                 os.path.exists(paired_counts_all) and 
                 os.path.exists(paired_counts_unique) and
+                os.path.exists(paired_rpkm_all) and
+                os.path.exists(paired_rpkm_unique) and
                 os.path.exists(all_bw) and 
                 os.path.exists(unique_bw)
             )
         else:  # single-end samples
             # Check if counts files exist
             single_counts_all = os.path.join(config['results_dir'], effective_name, f"{effective_name}_counts_single_all.txt")
             single_counts_unique = os.path.join(config['results_dir'], effective_name, f"{effective_name}_counts_single_unique.txt")
+
+            # Check if RPKM files exist
+            single_rpkm_all = os.path.join(config['results_dir'], effective_name, f"{effective_name}_rpkm_single_all.txt")
+            single_rpkm_unique = os.path.join(config['results_dir'], effective_name, f"{effective_name}_rpkm_single_unique.txt")
 
             # Check if coverage files exist
             all_bw = os.path.join(config['results_dir'], effective_name, f"{effective_name}_all.bw")
@@ -406,6 +416,8 @@ def is_sample_completed(config, sample, samples_df):
             final_results_exist = (
                 os.path.exists(single_counts_all) and 
                 os.path.exists(single_counts_unique) and
+                os.path.exists(single_rpkm_all) and
+                os.path.exists(single_rpkm_unique) and
                 os.path.exists(all_bw) and 
                 os.path.exists(unique_bw)
             )

rules/scripts/calc_rpkm.py

@@ -0,0 +1,53 @@
+import os
+
+import numpy as np
+import pandas as pd
+
+
+def rpkm(counts: np.ndarray, lengths: np.ndarray) -> np.ndarray:
+    """Calculate reads per kilobase transcript per million reads."""
+    total_reads = np.sum(counts, axis=0)
+    denominator = total_reads[np.newaxis, :] * lengths[:, np.newaxis]
+    return np.divide(1e9 * counts, denominator, out=np.zeros_like(counts, dtype=float), where=denominator != 0)
+
+
+def _log(message: str) -> None:
+    if hasattr(snakemake, "log") and snakemake.log:
+        os.makedirs(os.path.dirname(snakemake.log[0]), exist_ok=True)
+        with open(snakemake.log[0], "a", encoding="utf-8") as log_handle:
+            log_handle.write(f"{message}\n")
+
+
+def main() -> None:
+    input_path = snakemake.input[0]
+    output_path = snakemake.output[0]
+    os.makedirs(os.path.dirname(output_path), exist_ok=True)
+
+    _log(f"Loading counts from {input_path}")
+    counts_df = pd.read_csv(input_path, sep="\t", comment="#")
+
+    if "Length" not in counts_df.columns:
+        raise ValueError(f"Missing 'Length' column in {input_path}")
+
+    length_index = counts_df.columns.get_loc("Length")
+    sample_columns = list(counts_df.columns[length_index + 1 :])
+
+    if not sample_columns:
+        raise ValueError(f"No sample columns found in {input_path}")
+
+    counts_matrix = counts_df[sample_columns].to_numpy(dtype=float)
+    lengths = counts_df["Length"].to_numpy(dtype=float)
+
+    _log("Calculating RPKM values")
+    rpkm_matrix = rpkm(counts_matrix, lengths)
+
+    for col_index, column in enumerate(sample_columns):
+        counts_df[column] = rpkm_matrix[:, col_index].astype(float)
+
+    _log(f"Writing RPKM output to {output_path}")
+    counts_df.to_csv(output_path, sep="\t", index=False)
+    _log("RPKM calculation complete")
+
+
+if __name__ == "__main__":
+    main()