From 2a499175f9d2ceb5349318276be645ab3e9a5ee4 Mon Sep 17 00:00:00 2001 From: mtinti Date: Mon, 26 Jan 2026 13:18:31 +0000 Subject: [PATCH 1/5] Handle zero denominators in RPKM calculation --- README.md | 2 + Snakefile | 8 ++- envs/rpkm.yaml | 8 +++ rules/common.smk | 9 ++- rules/counting/rpkm.smk | 127 +++++++++++++++++++++++++++++++++++++ rules/results.smk | 8 ++- rules/sample_utils.py | 12 ++++ rules/scripts/calc_rpkm.py | 52 +++++++++++++++ 8 files changed, 221 insertions(+), 5 deletions(-) create mode 100644 envs/rpkm.yaml create mode 100644 rules/counting/rpkm.smk create mode 100644 rules/scripts/calc_rpkm.py diff --git a/README.md b/README.md index 6d74385..9107d24 100644 --- a/README.md +++ b/README.md @@ -38,6 +38,8 @@ Coverage Tracks (bamCoverage) ↓ Feature Counting (featureCounts) ↓ +RPKM Calculation + ↓ Results & Cleanup * Step skipped for nanopore data diff --git a/Snakefile b/Snakefile index 952c34c..31b9436 100644 --- a/Snakefile +++ b/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.") \ No newline at end of file + print("You can resume the pipeline with --keep-going to continue from where it left off.") diff --git a/envs/rpkm.yaml b/envs/rpkm.yaml new file mode 100644 index 0000000..058608e --- /dev/null +++ b/envs/rpkm.yaml @@ -0,0 +1,8 @@ +name: rpkm +channels: + - conda-forge + - defaults +dependencies: + - python=3.11 + - pandas + - numpy diff --git a/rules/common.smk b/rules/common.smk index c8c4aa0..e8ae45c 100644 --- a/rules/common.smk +++ b/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") \ No newline at end of file + return get_results_path(f"{run_tag}/sample_cleanup_complete.txt") diff --git a/rules/counting/rpkm.smk b/rules/counting/rpkm.smk new file mode 100644 index 0000000..0108ed0 --- /dev/null +++ b/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: + "rules/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: + "rules/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: + "rules/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: + "rules/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} + """ diff --git a/rules/results.smk b/rules/results.smk index ec2816c..2d5f79b 100644 --- a/rules/results.smk +++ b/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 - """ \ No newline at end of file + """ diff --git a/rules/sample_utils.py b/rules/sample_utils.py index 8a7d1b8..11d974a 100644 --- a/rules/sample_utils.py +++ b/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,6 +394,8 @@ 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) ) @@ -397,6 +403,10 @@ def is_sample_completed(config, sample, samples_df): # 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) ) diff --git a/rules/scripts/calc_rpkm.py b/rules/scripts/calc_rpkm.py new file mode 100644 index 0000000..320d5a7 --- /dev/null +++ b/rules/scripts/calc_rpkm.py @@ -0,0 +1,52 @@ +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) + + counts_df.loc[:, sample_columns] = rpkm_matrix + + _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() From 2315e9222b93d03a84352824debbf8e1036145d9 Mon Sep 17 00:00:00 2001 From: mtinti Date: Mon, 26 Jan 2026 13:28:29 +0000 Subject: [PATCH 2/5] Fix RPKM script path in Snakemake rules --- rules/counting/rpkm.smk | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/rules/counting/rpkm.smk b/rules/counting/rpkm.smk index 0108ed0..3667b03 100644 --- a/rules/counting/rpkm.smk +++ b/rules/counting/rpkm.smk @@ -41,7 +41,7 @@ rule rpkm_paired_all: conda: "../../envs/rpkm.yaml" script: - "rules/scripts/calc_rpkm.py" + "../scripts/calc_rpkm.py" rule rpkm_paired_unique: @@ -56,7 +56,7 @@ rule rpkm_paired_unique: conda: "../../envs/rpkm.yaml" script: - "rules/scripts/calc_rpkm.py" + "../scripts/calc_rpkm.py" rule rpkm_paired_complete: @@ -90,7 +90,7 @@ rule rpkm_single_all: conda: "../../envs/rpkm.yaml" script: - "rules/scripts/calc_rpkm.py" + "../scripts/calc_rpkm.py" rule rpkm_single_unique: @@ -105,7 +105,7 @@ rule rpkm_single_unique: conda: "../../envs/rpkm.yaml" script: - "rules/scripts/calc_rpkm.py" + "../scripts/calc_rpkm.py" rule rpkm_single_complete: From b8607d6b6d94c70de08ecac7a6d76c3b6d5825d6 Mon Sep 17 00:00:00 2001 From: mtinti Date: Mon, 26 Jan 2026 13:36:55 +0000 Subject: [PATCH 3/5] Cast count columns to float before RPKM assignment --- rules/scripts/calc_rpkm.py | 1 + 1 file changed, 1 insertion(+) diff --git a/rules/scripts/calc_rpkm.py b/rules/scripts/calc_rpkm.py index 320d5a7..a1dcdaf 100644 --- a/rules/scripts/calc_rpkm.py +++ b/rules/scripts/calc_rpkm.py @@ -35,6 +35,7 @@ def main() -> None: if not sample_columns: raise ValueError(f"No sample columns found in {input_path}") + counts_df.loc[:, sample_columns] = counts_df[sample_columns].astype(float) counts_matrix = counts_df[sample_columns].to_numpy(dtype=float) lengths = counts_df["Length"].to_numpy(dtype=float) From d3ff548e717997df5055be15a4ffd142e533b192 Mon Sep 17 00:00:00 2001 From: mtinti Date: Mon, 26 Jan 2026 13:46:13 +0000 Subject: [PATCH 4/5] Assign RPKM values per column to avoid pandas dtype errors --- rules/scripts/calc_rpkm.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rules/scripts/calc_rpkm.py b/rules/scripts/calc_rpkm.py index a1dcdaf..e40b2b7 100644 --- a/rules/scripts/calc_rpkm.py +++ b/rules/scripts/calc_rpkm.py @@ -35,14 +35,14 @@ def main() -> None: if not sample_columns: raise ValueError(f"No sample columns found in {input_path}") - counts_df.loc[:, sample_columns] = counts_df[sample_columns].astype(float) 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) - counts_df.loc[:, sample_columns] = rpkm_matrix + 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) From 273871022fe73c95d5bbc0a4420f79f6d799a526 Mon Sep 17 00:00:00 2001 From: mtinti Date: Mon, 26 Jan 2026 14:04:38 +0000 Subject: [PATCH 5/5] Show RPKM output in CI test log --- .github/workflows/test-PR.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/test-PR.yml b/.github/workflows/test-PR.yml index cfc9050..651ba3a 100644 --- a/.github/workflows/test-PR.yml +++ b/.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()