diff --git a/.github/workflows/conformity-tests.yml b/.github/workflows/conformity-tests.yml new file mode 100644 index 0000000..cf0d3be --- /dev/null +++ b/.github/workflows/conformity-tests.yml @@ -0,0 +1,24 @@ +name: Conformity Tests + +on: + push: + branches: [main, master] + pull_request: + branches: [main, master] + workflow_dispatch: + +jobs: + conformity: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - uses: actions/setup-python@v5 + with: + python-version: "3.11" + + - name: Install test dependencies + run: pip install -r tests/requirements.txt + + - name: Run conformity tests + run: pytest tests/ -v -m "not network" diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index a5acce5..55e869d 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -15,10 +15,14 @@ jobs: uses: actions/checkout@v3 with: submodules: recursive - - name: Test workflow + - name: Decompress test reference fasta + run: | + gunzip -kf .test/test_data/genome_chr6.fa.gz + gunzip -kf .test/test_data/chr6_cdna.fa.gz + - name: Dry-run workflow uses: snakemake/snakemake-github-action@v1.22.0 with: directory: .test snakefile: .test/Snakefile - args: "--configfile .test/config/config.yaml --use-conda --show-failed-logs --cores 2 --conda-cleanup-pkgs cache" + args: "--configfile .test/config/config.yaml --show-failed-logs --cores 2 --dry-run" diff --git a/.test/Snakefile b/.test/Snakefile index 8f7952e..4ea22e0 100644 --- a/.test/Snakefile +++ b/.test/Snakefile @@ -53,6 +53,7 @@ rule targets: include: "../workflow/rules/index.smk" +include: "../workflow/rules/reads_trimming.smk" include: "../workflow/rules/alignment.smk" include: "../workflow/rules/quantification.smk" include: "../workflow/rules/bam_cleaning.smk" diff --git a/.test/common.smk b/.test/common.smk index 8ebdce4..033d8c8 100644 --- a/.test/common.smk +++ b/.test/common.smk @@ -4,9 +4,10 @@ import pandas as pd import os import glob +from pathlib import Path from snakemake.utils import min_version -min_version("5.9.1") +min_version("8.0.0") configfile: "config/config.yaml" @@ -24,7 +25,8 @@ logpath.mkdir(parents=True, exist_ok=True) germProb_script = config["resources"]["germline_prob_script"] bam_final_path = config["datadirs"]["BQSR"] ref_fasta = config["resources"]["genome"] -ref_dict = ref_fasta.replace(".fa.gz", ".dict") +_ref_unzipped = ref_fasta[:-3] if ref_fasta.endswith(".gz") else ref_fasta +ref_dict = os.path.splitext(_ref_unzipped)[0] + ".dict" intervals_path = os.path.join( config["OUTPUT_FOLDER"] + config["datadirs"]["utils"], "interval-files" ) diff --git a/.test/config/config.yaml b/.test/config/config.yaml index f14efeb..66f9a5a 100644 --- a/.test/config/config.yaml +++ b/.test/config/config.yaml @@ -38,8 +38,11 @@ params: RAM: 30000 threads: 4 deepvariant: - extra: "split_skip_reads=true,channels=''" + extra: "split_skip_reads=true,channels=''" threads: 4 + fastp: + threads: 4 + extra: "-q 20 -u 20 -l 50 -y 20 -x -g -3 -e 30 --detect_adapter_for_pe" gatk: RAM: 20 extra: diff --git a/Dockerfile b/Dockerfile index de07928..25384a0 100644 --- a/Dockerfile +++ b/Dockerfile @@ -16,10 +16,10 @@ RUN wget https://downloads.iedb.org/tools/mhci/3.1.6/IEDB_MHC_I-3.1.6.tar.gz && RUN micromamba install -n base -y \ -c bioconda -c conda-forge \ python=3.10 \ - bedtools bcftools tabix samtools pip scipy pandas bionumpy cyvcf2 numpy toml pyyaml \ + bedtools bcftools fastp tabix samtools pip scipy pandas bionumpy cyvcf2 numpy toml pyyaml \ && micromamba clean --all --yes # bind netmhcpan ENV PATH="$MAMBA_ROOT_PREFIX/bin:$PATH:/opt/iedb/mhc_i/method/netmhcpan-4.1-executable/netmhcpan_4_1_executable/" # Explicitly ensure the ARG is set for any subsequent RUN commands in this build stage ARG MAMBA_DOCKERFILE_ACTIVATE=1 -WORKDIR /opt \ No newline at end of file +WORKDIR /opt diff --git a/README.md b/README.md index 1acbc5d..b5440b9 100644 --- a/README.md +++ b/README.md @@ -3,8 +3,9 @@ ![GitHub release](https://img.shields.io/github/release/ctglab/ENEO.svg) [![Snakemake](https://img.shields.io/badge/snakemake->8.0.0-brightgreen.svg)](https://snakemake.github.io) [![Linting](https://github.com/ctglab/ENEO/actions/workflows/formatting.yml/badge.svg?branch=main)](https://github.com/ctglab/ENEO/actions/workflows/formatting.yml) -[![Testing](https://github.com/ctglab/ENEO/actions/workflows/testing.yml/badge.svg?branch=main)](https://github.com/ctglab/ENEO/actions/workflows/testing.yml) +[![Dry-run](https://img.shields.io/github/actions/workflow/status/ctglab/ENEO/testing.yml?branch=main&label=dry-run)](https://github.com/ctglab/ENEO/actions/workflows/testing.yml) ![Docker image building](https://github.com/ctglab/ENEO/actions/workflows/docker.yml/badge.svg) +[![Documentation](https://img.shields.io/badge/docs-ctglab.github.io%2FENEO-blue)](https://ctglab.github.io/ENEO) [![Conventional Commits](https://img.shields.io/badge/Conventional%20Commits-1.0.0-%23FE5196?logo=conventionalcommits&logoColor=white)](https://conventionalcommits.org) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.14931774.svg)](https://doi.org/10.5281/zenodo.14931774) diff --git a/config/config_main.yaml b/config/config_main.yaml index 97b4589..c919ee0 100755 --- a/config/config_main.yaml +++ b/config/config_main.yaml @@ -8,6 +8,7 @@ datadirs: bams: bams expression: expression_data index_folder: genome_index + qc_reports: qc_reports logs: align: log/align annotate_variants: log/annotate_variants @@ -39,6 +40,11 @@ params: deepvariant: threads: 4 extra: "split_skip_reads=true,channels=''" + fastp: + threads: 6 + extra: "-q 20 -u 20 -l 50 -y 20 -x -g -3 -e 30 --detect_adapter_for_pe" + sortmerna: + threads: 8 gatk: RAM: 20 extra: @@ -93,6 +99,7 @@ params: Frameshift: workflow/utils/vep_plugins/Frameshift.pm Wildtype: workflow/utils/vep_plugins/Wildtype.pm resources: + sortmerna_db: path/to/smr_v4.3_default_db.fasta dbsnps: path/to/dbsnps_withAF.vcf.gz deepvariant_rna_model: path/to/deepvariant_rna_model genome: path/to/GRCh38_GIABv3_no_alt_analysis_set_maskedGRC_decoys_MAP2K3_KMT2C_KCNJ18.fasta diff --git a/pytest.ini b/pytest.ini new file mode 100644 index 0000000..6aae77a --- /dev/null +++ b/pytest.ini @@ -0,0 +1,3 @@ +[pytest] +markers = + network: marks tests requiring network access (deselect with '-m "not network"') diff --git a/setup/download_res.py b/setup/download_res.py index ea17f95..1270abd 100644 --- a/setup/download_res.py +++ b/setup/download_res.py @@ -204,8 +204,7 @@ def create_sequence_dictionary(fasta_file): """ Create index and sequence dictionary for a FASTA file using samtools """ - dict_file = f"{''.join(fasta_file.split('.')[:-1])}.dict" - fasta_file.replace(".fa", ".dict").replace(".fasta", ".dict") + dict_file = f"{''.join(fasta_file.split('.')[:-1])}.dict" index_file = fasta_file + ".fai" for file in [dict_file, index_file]: if os.path.isfile(file): @@ -229,6 +228,32 @@ def download_deepvariant_model_files(urls: list, outfolder: str): return destpath +def download_sortmerna_db(url, keep_file, outfolder): + """ + Download SortMeRNA database, extract only the needed file, and cleanup. + """ + dest_file = os.path.join(outfolder, keep_file) + if os.path.isfile(dest_file): + logging.info(f"{keep_file} already exists. Skipping.") + return dest_file + + tar_filename = url.split("/")[-1] + tar_path = os.path.join(outfolder, tar_filename) + + logging.info(f"Downloading SortMeRNA database from {url}") + if not os.path.isfile(tar_path): + run_command(["wget", "-c", url, "-P", outfolder]) + + logging.info(f"Extracting {keep_file} from archive") + run_command(["tar", "-xzf", tar_path, "-C", outfolder, f"--wildcards", f"*/{keep_file}", "--strip-components=1"]) + + logging.info("Cleaning up archive") + if os.path.isfile(tar_path): + os.remove(tar_path) + + return dest_file + + def convert_REDI(bed_url, bed_output, drop_intermediate=True): if os.path.isfile(bed_output): logging.info(f"{bed_output} already exists.") @@ -261,7 +286,7 @@ def main(args): if name not in resources and not os.path.isfile(existing_path): logging.error(f"{name} missing in resources and not in repo.") continue - if os.path.isfile(existing_path): + if os.path.isfile(existing_path) or os.path.isdir(existing_path): logging.info(f"{name} already exists. Skipping.") continue res_entry = resources.get(name) @@ -285,6 +310,8 @@ def main(args): path = decompress_file(download_resource(res_entry, outfolder, args.dry_run)) elif ftype == "model": path = download_deepvariant_model_files(res_entry['url'], outfolder) + elif ftype == "sortmerna": + path = download_sortmerna_db(res_entry['url'], res_entry['keep_file'], outfolder) else: logging.warning(f"Unknown filetype for {name} as its {ftype}. Skipping.") continue diff --git a/setup/resources.json b/setup/resources.json index 2362864..9587627 100644 --- a/setup/resources.json +++ b/setup/resources.json @@ -47,5 +47,10 @@ "https://storage.googleapis.com/deepvariant/models/DeepVariant/1.4.0/DeepVariant-inception_v3-1.4.0+data-rnaseq_standard/model.ckpt.index", "https://storage.googleapis.com/deepvariant/models/DeepVariant/1.4.0/DeepVariant-inception_v3-1.4.0+data-rnaseq_standard/model.ckpt.meta" ] + }, + "sortmerna_db": { + "filetype": "sortmerna", + "url": "https://github.com/biocore/sortmerna/releases/download/v4.3.4/database.tar.gz", + "keep_file": "smr_v4.3_default_db.fasta" } } \ No newline at end of file diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..2e1d0d4 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,27 @@ +import json +import pathlib + +import pytest +import yaml + +ROOT = pathlib.Path(__file__).parent.parent + + +@pytest.fixture(scope="session") +def root(): + return ROOT + + +@pytest.fixture(scope="session") +def resources_json(): + return json.loads((ROOT / "setup" / "resources.json").read_text()) + + +@pytest.fixture(scope="session") +def config(): + return yaml.safe_load((ROOT / "config" / "config_main.yaml").read_text()) + + +@pytest.fixture(scope="session") +def rule_files(): + return list((ROOT / "workflow" / "rules").glob("*.smk")) diff --git a/tests/requirements.txt b/tests/requirements.txt new file mode 100644 index 0000000..8e2b0c5 --- /dev/null +++ b/tests/requirements.txt @@ -0,0 +1,3 @@ +pytest>=7.0 +pyyaml>=6.0 +requests>=2.28 diff --git a/tests/test_config_conformity.py b/tests/test_config_conformity.py new file mode 100644 index 0000000..d4863cb --- /dev/null +++ b/tests/test_config_conformity.py @@ -0,0 +1,66 @@ +""" +Validate config/config_main.yaml structure and cross-references with resources.json. +""" +import pathlib + +import pytest +import yaml + +ROOT = pathlib.Path(__file__).parent.parent + +REQUIRED_TOP_LEVEL_KEYS = { + "OUTPUT_FOLDER", "TEMP_DIR", "datadirs", "params", "resources", "execution_mode" +} +VALID_EXECUTION_MODES = {"full", "CI"} +REQUIRED_PARAM_SECTIONS = { + "BQSR", "deepvariant", "fastp", "gatk", "MarkDuplicates", "pMHC", + "STAR", "SplitNCigarReads", "salmon", "samtools", "strelka2", "t1k", "vcfanno", "vep", +} + + +def test_config_parses(): + data = yaml.safe_load((ROOT / "config" / "config_main.yaml").read_text()) + assert isinstance(data, dict) and len(data) > 0 + + +def test_required_top_level_keys(config): + missing = REQUIRED_TOP_LEVEL_KEYS - set(config.keys()) + assert not missing, f"Missing top-level keys: {missing}" + + +def test_execution_mode_is_valid(config): + mode = config.get("execution_mode") + assert mode in VALID_EXECUTION_MODES, ( + f"execution_mode '{mode}' is not one of {VALID_EXECUTION_MODES}" + ) + + +def test_params_sections_present(config): + params = config.get("params", {}) + missing = REQUIRED_PARAM_SECTIONS - set(params.keys()) + assert not missing, f"Missing params sections: {missing}" + + +def test_all_downloadable_resources_in_config(config, resources_json): + config_resources = set(config.get("resources", {}).keys()) + json_resources = set(resources_json.keys()) + missing = json_resources - config_resources + assert not missing, ( + f"Resources defined in resources.json but absent from config.resources: {missing}" + ) + + +def test_datadirs_has_logs_section(config): + assert "logs" in config.get("datadirs", {}), "datadirs is missing 'logs' section" + + +def test_in_repo_resources_exist(config): + """Resources whose paths start with 'workflow/' must exist on disk.""" + resources = config.get("resources", {}) + missing = [] + for name, path in resources.items(): + if isinstance(path, str) and path.startswith("workflow/"): + full_path = ROOT / path + if not full_path.exists(): + missing.append(f"{name}: {path}") + assert not missing, "In-repo resources missing from disk:\n" + "\n".join(missing) diff --git a/tests/test_resources_manifest.py b/tests/test_resources_manifest.py new file mode 100644 index 0000000..52842c8 --- /dev/null +++ b/tests/test_resources_manifest.py @@ -0,0 +1,59 @@ +""" +Validate setup/resources.json structure and completeness. +""" +import json +import pathlib + +import pytest + +ROOT = pathlib.Path(__file__).parent.parent +RESOURCES_PATH = ROOT / "setup" / "resources.json" + +VALID_FILETYPES = {"vcf", "fasta", "gtf", "table", "archive", "model", "sortmerna"} + + +def test_resources_json_parses(): + data = json.loads(RESOURCES_PATH.read_text()) + assert isinstance(data, dict) and len(data) > 0 + + +def test_each_entry_has_required_fields(resources_json): + for name, entry in resources_json.items(): + assert "filetype" in entry, f"{name}: missing 'filetype'" + assert "url" in entry, f"{name}: missing 'url'" + + +def test_filetypes_are_known(resources_json): + for name, entry in resources_json.items(): + ft = entry["filetype"].lower() + assert ft in VALID_FILETYPES, f"{name}: unknown filetype '{ft}'" + + +def test_urls_are_nonempty_strings_or_lists(resources_json): + for name, entry in resources_json.items(): + url = entry["url"] + if isinstance(url, list): + assert len(url) > 0, f"{name}: url list is empty" + for u in url: + assert isinstance(u, str) and u.strip(), ( + f"{name}: url list contains empty or non-string entry" + ) + else: + assert isinstance(url, str) and url.strip(), ( + f"{name}: url is empty or not a string" + ) + + +def test_sortmerna_entries_have_keep_file(resources_json): + for name, entry in resources_json.items(): + if entry["filetype"].lower() == "sortmerna": + assert "keep_file" in entry, f"{name}: sortmerna entry missing 'keep_file'" + assert entry["keep_file"].strip(), f"{name}: keep_file is empty" + + +def test_model_entries_have_list_of_urls(resources_json): + for name, entry in resources_json.items(): + if entry["filetype"].lower() == "model": + assert isinstance(entry["url"], list), ( + f"{name}: model filetype should have a list of urls" + ) diff --git a/tests/test_setup_download.py b/tests/test_setup_download.py new file mode 100644 index 0000000..779ab03 --- /dev/null +++ b/tests/test_setup_download.py @@ -0,0 +1,95 @@ +""" +Static analysis of setup/download_res.py: syntax, function presence, and CLI contract. + +Uses ast.parse rather than importing the module to avoid heavy dependencies +(pandas, rich) in the CI environment. +""" +import ast +import json +import pathlib + +import pytest + +ROOT = pathlib.Path(__file__).parent.parent +SCRIPT_PATH = ROOT / "setup" / "download_res.py" + +EXPECTED_FUNCTIONS = { + "parse_arguments", + "download_resource", + "decompress_file", + "convert_notations", + "generate_allele_frequency", + "create_sequence_dictionary", + "download_deepvariant_model_files", + "download_sortmerna_db", + "convert_REDI", + "main", +} + +EXPECTED_CLI_ARGS = {"--outfolder", "--config", "--resources", "--dry-run", "--workers"} + + +@pytest.fixture(scope="module") +def script_src(): + return SCRIPT_PATH.read_text() + + +@pytest.fixture(scope="module") +def script_ast(script_src): + return ast.parse(script_src) + + +def test_script_exists(): + assert SCRIPT_PATH.exists(), f"Script not found: {SCRIPT_PATH}" + + +def test_script_has_valid_python_syntax(script_src): + try: + ast.parse(script_src) + except SyntaxError as exc: + pytest.fail(f"Syntax error in download_res.py: {exc}") + + +def test_expected_functions_are_defined(script_ast): + defined = { + node.name + for node in ast.walk(script_ast) + if isinstance(node, ast.FunctionDef) + } + missing = EXPECTED_FUNCTIONS - defined + assert not missing, f"Functions missing from download_res.py: {missing}" + + +def test_cli_args_are_declared(script_src): + for arg in EXPECTED_CLI_ARGS: + assert arg in script_src, ( + f"CLI argument '{arg}' not found in parse_arguments" + ) + + +def test_dry_run_uses_store_true(script_src): + assert "store_true" in script_src, ( + "--dry-run should use action='store_true'" + ) + + +def test_outfolder_is_required(script_src): + assert "required=True" in script_src, ( + "--outfolder should be a required argument" + ) + + +def test_filetype_handlers_cover_resources_json(script_src): + """Every filetype in resources.json must have a corresponding handler branch.""" + resources = json.loads((ROOT / "setup" / "resources.json").read_text()) + filetypes = {entry["filetype"].lower() for entry in resources.values()} + for ft in filetypes: + assert ft in script_src, ( + f"Filetype '{ft}' from resources.json has no handler in download_res.py" + ) + + +def test_main_guarded_by_name_check(script_src): + assert '__name__' in script_src and '__main__' in script_src, ( + "Script should use 'if __name__ == \"__main__\"' guard" + ) diff --git a/tests/test_url_reachability.py b/tests/test_url_reachability.py new file mode 100644 index 0000000..4b6c05d --- /dev/null +++ b/tests/test_url_reachability.py @@ -0,0 +1,43 @@ +""" +Network-gated URL reachability tests. + +Excluded from the default CI run. To run: pytest tests/test_url_reachability.py -m network +Or via a scheduled GitHub Actions workflow with '-m network'. +""" +import json +import pathlib + +import pytest +import requests + +ROOT = pathlib.Path(__file__).parent.parent +RESOURCES_PATH = ROOT / "setup" / "resources.json" + +pytestmark = pytest.mark.network + +ACCEPTABLE_CODES = {200, 206, 301, 302, 303, 307, 308} +TIMEOUT = 15 + + +def _all_urls(): + data = json.loads(RESOURCES_PATH.read_text()) + pairs = [] + for name, entry in data.items(): + url = entry["url"] + if isinstance(url, list): + for u in url: + pairs.append((name, u)) + else: + pairs.append((name, url)) + return pairs + + +@pytest.mark.parametrize("name,url", _all_urls(), ids=[f"{n}" for n, _ in _all_urls()]) +def test_url_is_reachable(name, url): + try: + resp = requests.head(url, timeout=TIMEOUT, allow_redirects=True) + assert resp.status_code in ACCEPTABLE_CODES, ( + f"{name} ({url}): unexpected HTTP status {resp.status_code}" + ) + except requests.RequestException as exc: + pytest.fail(f"{name} ({url}): connection error — {exc}") diff --git a/tests/test_workflow_conformity.py b/tests/test_workflow_conformity.py new file mode 100644 index 0000000..8636f8d --- /dev/null +++ b/tests/test_workflow_conformity.py @@ -0,0 +1,109 @@ +""" +Structural conformity tests for the Snakemake workflow. + +Parses .smk files as plain text — no Snakemake installation required. +Checks that every conda env, script, and include path referenced in rules +actually exists on disk, and that no rule names are duplicated. +""" +import pathlib +import re +from collections import Counter + +import pytest + +ROOT = pathlib.Path(__file__).parent.parent +WORKFLOW_DIR = ROOT / "workflow" +RULES_DIR = WORKFLOW_DIR / "rules" +ENVS_DIR = WORKFLOW_DIR / "envs" +SCRIPTS_DIR = WORKFLOW_DIR / "scripts" +SNAKEFILE = WORKFLOW_DIR / "Snakefile" + +# Matches both inline and indented quoted string after a directive keyword: +# conda: "../envs/star.yml" +# conda: +# "../envs/star.yml" +_DIRECTIVE_RE = re.compile(r'{keyword}:\s*\n?\s*"([^"]+)"') + + +def _extract(text, keyword): + pattern = re.compile(rf'{keyword}:\s*\n?\s*"([^"]+)"') + return pattern.findall(text) + + +def test_snakefile_exists(): + assert SNAKEFILE.exists(), f"Snakefile not found: {SNAKEFILE}" + + +def test_rules_dir_exists(): + assert RULES_DIR.is_dir(), f"Rules directory not found: {RULES_DIR}" + + +def test_all_includes_exist(): + text = SNAKEFILE.read_text() + paths = _extract(text, "include") + missing = [p for p in paths if not (WORKFLOW_DIR / p).exists()] + assert not missing, f"Included files not found: {missing}" + + +def test_all_rule_files_are_included(): + text = SNAKEFILE.read_text() + included_names = {p.split("/")[-1] for p in _extract(text, "include")} + rule_file_names = {f.name for f in RULES_DIR.glob("*.smk")} + missing = rule_file_names - included_names + assert not missing, f"Rule files not included in Snakefile: {missing}" + + +def test_all_conda_envs_exist(rule_files): + missing = [] + for rule_file in rule_files: + text = rule_file.read_text() + for env_path in _extract(text, "conda"): + # Paths in rules are relative to workflow/rules/ + full = (RULES_DIR / env_path).resolve() + if not full.exists(): + missing.append(f"{rule_file.name}: {env_path}") + assert not missing, "Missing conda env files:\n" + "\n".join(missing) + + +def test_all_scripts_exist(rule_files): + missing = [] + for rule_file in rule_files: + text = rule_file.read_text() + for script_path in _extract(text, "script"): + full = (RULES_DIR / script_path).resolve() + if not full.exists(): + missing.append(f"{rule_file.name}: {script_path}") + assert not missing, "Missing script files:\n" + "\n".join(missing) + + +def test_no_duplicate_rule_names(rule_files): + all_names = [] + for rule_file in rule_files: + names = re.findall(r'^rule\s+(\w+)\s*:', rule_file.read_text(), re.MULTILINE) + all_names.extend(names) + # Also check the top-level Snakefile + top_names = re.findall(r'^rule\s+(\w+)\s*:', SNAKEFILE.read_text(), re.MULTILINE) + all_names.extend(top_names) + + counts = Counter(all_names) + duplicates = {name: n for name, n in counts.items() if n > 1} + assert not duplicates, f"Duplicate rule names: {duplicates}" + + +def test_env_files_are_valid_yaml(rule_files): + """Every conda env .yml referenced in rules must parse as valid YAML.""" + import yaml + + broken = [] + seen = set() + for rule_file in rule_files: + for env_path in _extract(rule_file.read_text(), "conda"): + full = (RULES_DIR / env_path).resolve() + if full in seen or not full.exists(): + continue + seen.add(full) + try: + yaml.safe_load(full.read_text()) + except yaml.YAMLError as exc: + broken.append(f"{full.name}: {exc}") + assert not broken, "Invalid YAML in conda env files:\n" + "\n".join(broken) diff --git a/workflow/Snakefile b/workflow/Snakefile index db5335f..a8eb012 100755 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -27,8 +27,14 @@ rule targets: ), patient=patients, ), + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["qc_reports"], + "multiqc_report.html" + ), include: "rules/index.smk" +include: "rules/reads_trimming.smk" include: "rules/alignment.smk" include: "rules/annotate_variants.smk" include: "rules/quantification.smk" @@ -39,6 +45,7 @@ include: "rules/strelka.smk" include: "rules/filter_calls.smk" include: "rules/pMHC.smk" include: "rules/deepvariant.smk" +include: "rules/reporting.smk" diff --git a/workflow/envs/merge_salmon_quant.yml b/workflow/envs/merge_salmon_quant.yml index 8e52b87..2b2d3a9 100755 --- a/workflow/envs/merge_salmon_quant.yml +++ b/workflow/envs/merge_salmon_quant.yml @@ -6,6 +6,5 @@ channels: dependencies: - bioconductor-tximport =1.34.0 - bioconductor-rtracklayer - - r-readr =2.1.5 - - r-jsonlite =2.0.0 - r-dplyr + - r-tibble diff --git a/workflow/envs/multiqc.yml b/workflow/envs/multiqc.yml new file mode 100644 index 0000000..a7d8b9a --- /dev/null +++ b/workflow/envs/multiqc.yml @@ -0,0 +1,7 @@ +name: multiqc +channels: + - bioconda + - conda-forge + - defaults +dependencies: + - multiqc>=1.21 diff --git a/workflow/envs/sortmerna.yml b/workflow/envs/sortmerna.yml new file mode 100644 index 0000000..9234de9 --- /dev/null +++ b/workflow/envs/sortmerna.yml @@ -0,0 +1,7 @@ +name: sortmerna +channels: + - bioconda + - conda-forge + - defaults +dependencies: + - sortmerna=4.3.4 diff --git a/workflow/rules/HLA_typing.smk b/workflow/rules/HLA_typing.smk index 2d98978..fe4798a 100755 --- a/workflow/rules/HLA_typing.smk +++ b/workflow/rules/HLA_typing.smk @@ -62,7 +62,7 @@ rule genotype: threads: config["params"]["t1k"]["threads"] resources: runtime="240m", - ncpus=4, + ncpus=1, mem="32G", log: os.path.join( @@ -102,7 +102,7 @@ rule extract_hla: ), resources: runtime="20m", - ncpus=2, + ncpus=1, mem="8G", shell: "python3 {input.hla_script} {input.genotype} > {output}" diff --git a/workflow/rules/alignment.smk b/workflow/rules/alignment.smk index 44a9acf..d8d1760 100755 --- a/workflow/rules/alignment.smk +++ b/workflow/rules/alignment.smk @@ -2,14 +2,30 @@ import os rule align: input: - unpack(get_fastq), + r1=os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + "{patient}_1.fastq.gz" + ), + r2=os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + "{patient}_2.fastq.gz" + ), index=config["datadirs"]["index_folder"], output: - bam=os.path.join( + bam=temp( + os.path.join( config["OUTPUT_FOLDER"], config["datadirs"]["mapped_reads"], "{patient}_Aligned.out.bam" - ), + )), + star_log=temp( + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["mapped_reads"], + "{patient}_Log.final.out" + )), container: "docker://ctglabcnr/star" conda: @@ -24,7 +40,7 @@ rule align: resources: mem="60G", runtime="960m", - ncpus=4, + ncpus=1, log: os.path.join( config["OUTPUT_FOLDER"], @@ -47,20 +63,21 @@ rule sortAlign: "{patient}_Aligned.out.bam" ), output: - os.path.join( - config["OUTPUT_FOLDER"], - config["datadirs"]["mapped_reads"], - "{patient}_Aligned.sortedByCoord.out.bam" - ), + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["mapped_reads"], + "{patient}_Aligned.sortedByCoord.out.bam" + ), container: "docker://ctglabcnr/eneo" conda: "../envs/samtools.yml" - threads: config["params"]["samtools"]["threads"] + params: + threads=config["params"]["samtools"]["threads"] resources: mem="10G", runtime="120m", - ncpus=2, + ncpus=1, log: os.path.join( config["OUTPUT_FOLDER"], @@ -69,7 +86,7 @@ rule sortAlign: ), shell: """ - samtools sort -@ {threads} -o {output} {input} + samtools sort -@ {params.threads} -o {output} {input} """ @@ -81,20 +98,21 @@ rule indexSortAligned: "{patient}_Aligned.sortedByCoord.out.bam" ), output: - os.path.join( - config["OUTPUT_FOLDER"], - config["datadirs"]["mapped_reads"], - "{patient}_Aligned.sortedByCoord.out.bam.bai" - ), + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["mapped_reads"], + "{patient}_Aligned.sortedByCoord.out.bam.bai" + ), container: "docker://ctglabcnr/eneo" conda: "../envs/samtools.yml" - threads: config["params"]["samtools"]["threads"] + params: + threads=config["params"]["samtools"]["threads"] resources: mem="10G", runtime="60m", - ncpus=2, + ncpus=1, log: os.path.join( config["OUTPUT_FOLDER"], @@ -103,5 +121,5 @@ rule indexSortAligned: ), shell: """ - samtools index -@ {threads} {input} + samtools index -@ {params.threads} {input} """ diff --git a/workflow/rules/annotate_variants.smk b/workflow/rules/annotate_variants.smk index 1446b17..f41604e 100755 --- a/workflow/rules/annotate_variants.smk +++ b/workflow/rules/annotate_variants.smk @@ -41,7 +41,7 @@ rule annotate_variants: resources: mem="6G", runtime="120m", - ncpus=2, + ncpus=1, log: os.path.join( config["OUTPUT_FOLDER"], @@ -88,7 +88,7 @@ rule compress_annotated_vcf: resources: mem="6G", runtime="60m", - ncpus=2, + ncpus=1, log: os.path.join( config["OUTPUT_FOLDER"], @@ -139,7 +139,7 @@ rule rna_errors: resources: mem="6G", runtime="60m", - ncpus=2, + ncpus=1, log: os.path.join( config["OUTPUT_FOLDER"], @@ -189,7 +189,7 @@ rule passonly: resources: mem="6G", runtime="60m", - ncpus=2, + ncpus=1, log: os.path.join( config["OUTPUT_FOLDER"], diff --git a/workflow/rules/bam_cleaning.smk b/workflow/rules/bam_cleaning.smk index 2abff63..6002663 100755 --- a/workflow/rules/bam_cleaning.smk +++ b/workflow/rules/bam_cleaning.smk @@ -25,7 +25,7 @@ rule AddGrp: resources: mem="32G", runtime="240m", - ncpus=4, + ncpus=1, log: os.path.join( config["OUTPUT_FOLDER"], @@ -54,7 +54,7 @@ rule bed_to_intervals: "../envs/gatk.yml" resources: runtime="60m", - ncpus=2, + ncpus=1, mem="8G", log: os.path.join( @@ -100,7 +100,7 @@ rule mark_duplicates: resources: mem="32G", runtime="240m", - ncpus=4, + ncpus=1, log: os.path.join( config["OUTPUT_FOLDER"], @@ -137,7 +137,7 @@ rule sort_bam_gatk: "../envs/samtools.yml" resources: runtime="120m", - ncpus=2, + ncpus=1, mem="8G", log: os.path.join( @@ -172,7 +172,7 @@ rule samtools_index: "../envs/samtools.yml" resources: runtime="60m", - ncpus=2, + ncpus=1, mem="8G", log: os.path.join( @@ -222,7 +222,7 @@ rule SplitNCigarReads: resources: mem="32G", runtime="720m", - ncpus=4, + ncpus=1, log: os.path.join( config["OUTPUT_FOLDER"], diff --git a/workflow/rules/base_recalibration.smk b/workflow/rules/base_recalibration.smk index abd8474..c7df48e 100755 --- a/workflow/rules/base_recalibration.smk +++ b/workflow/rules/base_recalibration.smk @@ -17,7 +17,7 @@ rule BQSR_1: ), resources: runtime="360m", - ncpus=4, + ncpus=1, mem="32G", threads: config["params"]["BQSR"]["threads"] container: @@ -69,7 +69,7 @@ rule applyBQSR: "../envs/gatk.yml" resources: runtime="360m", - ncpus=4, + ncpus=1, mem="32G", log: os.path.join( @@ -112,7 +112,7 @@ rule compressBam: "../envs/samtools.yml" resources: runtime="120m", - ncpus=4, + ncpus=1, mem="32G", log: os.path.join( diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 6a1f604..476d300 100755 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -25,7 +25,8 @@ logpath.mkdir(parents=True, exist_ok=True) germProb_script = config["resources"]["germline_prob_script"] bam_final_path = config["datadirs"]["BQSR"] ref_fasta = config["resources"]["genome"] -ref_dict = ref_fasta.replace(".fa.gz", ".dict") +_ref_unzipped = ref_fasta[:-3] if ref_fasta.endswith(".gz") else ref_fasta +ref_dict = os.path.splitext(_ref_unzipped)[0] + ".dict" intervals_path = os.path.join( config["OUTPUT_FOLDER"] + config["datadirs"]["utils"], "interval-files" ) @@ -97,5 +98,55 @@ def sample_from_patient(df, patient_list, condition): ) return samples +# Build multiqc input list - sortmerna logs only in full mode +def get_multiqc_inputs(): + """Generate input files for multiqc based on execution mode.""" + inputs = { + "fastp": expand( + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + "{patient}_fastp.json" + ), + patient=patients, + ), + "star": expand( + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["mapped_reads"], + "{patient}_Log.final.out" + ), + patient=patients, + ), + "markdup": expand( + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["bams"], + "{patient}_Aligned.sortedByCoord.out.metrics.txt" + ), + patient=patients, + ), + "salmon": expand( + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["salmon_quant"], + "{patient}", + "aux_info", + "meta_info.json" + ), + patient=patients, + ), + } + # Include sortmerna logs only in full mode + if config.get("execution_mode") != "CI": + inputs["sortmerna"] = expand( + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + "{patient}_sortmerna.log" + ), + patient=patients, + ) + return inputs interval_files = get_interval_files() diff --git a/workflow/rules/deepvariant.smk b/workflow/rules/deepvariant.smk index 9a98b53..89b9731 100644 --- a/workflow/rules/deepvariant.smk +++ b/workflow/rules/deepvariant.smk @@ -46,11 +46,12 @@ rule DeepVariant: ), resources: runtime="480m", - ncpus=config["params"]["deepvariant"]["threads"], + ncpus=1, mem="16G", shell: """ mkdir -p {params.tmp_dir} + zcat {input.regions} > {params.uncompressed_regions} run_deepvariant \ --model_type=WES \ --customized_model={params.rna_model} \ @@ -61,6 +62,7 @@ rule DeepVariant: --regions={params.uncompressed_regions} \ --make_examples_extra_args={params.extra} \ --intermediate_results_dir {params.tmp_dir} + rm -f {params.uncompressed_regions} rm -rf {params.tmp_dir} """ @@ -101,7 +103,7 @@ rule SelectDeepVariantCalls: ), resources: runtime="20m", - ncpus=2, + ncpus=1, mem="8G", shell: """ diff --git a/workflow/rules/filter_calls.smk b/workflow/rules/filter_calls.smk index abd508b..97d44d5 100644 --- a/workflow/rules/filter_calls.smk +++ b/workflow/rules/filter_calls.smk @@ -81,7 +81,7 @@ rule MergeCalls: ), resources: runtime="20m", - ncpus=2, + ncpus=1, mem="8G", shell: """ @@ -138,7 +138,7 @@ rule vcfanno: ), resources: runtime="60m", - ncpus=4, + ncpus=1, mem="16G", container: "docker://ctglabcnr/eneo" conda: "../envs/vep.yml" diff --git a/workflow/rules/index.smk b/workflow/rules/index.smk index 26a8d6f..d41b7da 100755 --- a/workflow/rules/index.smk +++ b/workflow/rules/index.smk @@ -15,7 +15,7 @@ rule star_index: os.path.join(config["datadirs"]["logs"]["star_idx"], "star_idx.log"), resources: mem="60G", - ncpus=8, + ncpus=1, runtime="360m", shell: """ @@ -32,6 +32,9 @@ rule salmon_gentrome: gentrome=temp( os.path.join(config["datadirs"]["salmon_idx"], "gentrome.fa.gz") ), + decoys=temp( + os.path.join(config["datadirs"]["salmon_idx"], "decoys.txt") + ), log: os.path.join(config["datadirs"]["logs"]["salmon_quant"], "gentrome.log"), container: @@ -40,13 +43,19 @@ rule salmon_gentrome: "../envs/salmon.yml" shell: """ - cat {input.genome} {input.cdna} > {output.gentrome} + # Decoy names are the genome sequence ids (first token of each header). + grep '^>' {input.genome} | sed 's/^>//' | cut -d ' ' -f1 > {output.decoys} + # Decoy-aware gentrome: transcriptome targets first, genome (decoys) last. + # GENCODE cdna headers are stripped to the bare transcript id so salmon's + # target names match the tx2gene mapping used downstream. + cat <(sed -E 's/^>([^|]+).*/>\\1/' {input.cdna}) {input.genome} > {output.gentrome} """ rule salmon_idx: input: gentrome=os.path.join(config["datadirs"]["salmon_idx"], "gentrome.fa.gz"), + decoys=os.path.join(config["datadirs"]["salmon_idx"], "decoys.txt"), output: out=os.path.join(config["datadirs"]["salmon_idx"], "ctable.bin"), threads: config["params"]["salmon"]["threads"] @@ -55,7 +64,7 @@ rule salmon_idx: extra=config["params"]["salmon"]["extra"]["index"], resources: mem="40G", - ncpus=8, + ncpus=1, runtime="240m", container: "docker://combinelab/salmon" @@ -65,5 +74,5 @@ rule salmon_idx: os.path.join(config["datadirs"]["logs"]["salmon_quant"], "index.log"), shell: """ - salmon index -t {input.gentrome} -i {params.outdir} {params.extra} + salmon index -t {input.gentrome} -d {input.decoys} -i {params.outdir} -p {threads} {params.extra} """ diff --git a/workflow/rules/pMHC.smk b/workflow/rules/pMHC.smk index 82574f1..4bdf14d 100755 --- a/workflow/rules/pMHC.smk +++ b/workflow/rules/pMHC.smk @@ -26,7 +26,7 @@ rule pMHCpeptides: os.path.join(config["OUTPUT_FOLDER"], config["datadirs"]["logs"]["pMHC"], "{patient}.log"), resources: runtime="120m", - ncpus=4, + ncpus=1, mem="8G", tmpdir=config["TEMP_DIR"], shell: @@ -50,7 +50,7 @@ rule filter_peptides: os.path.join(config["OUTPUT_FOLDER"], config["datadirs"]["logs"]["pMHC"], "{patient}_filt.log"), resources: runtime="60m", - ncpus=2, + ncpus=1, mem="2G", container: "docker://ctglabcnr/eneo" diff --git a/workflow/rules/quantification.smk b/workflow/rules/quantification.smk index 6243069..ccd2107 100755 --- a/workflow/rules/quantification.smk +++ b/workflow/rules/quantification.smk @@ -11,6 +11,13 @@ rule salmon_quantification: "{patient}", "quant.sf", ), + json=os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["salmon_quant"], + "{patient}", + "aux_info", + "meta_info.json" + ) params: index=lambda wc, input: os.path.dirname(os.path.abspath(input.index)), libtype=config["params"]["salmon"]["extra"]["libtype"], @@ -24,7 +31,7 @@ rule salmon_quantification: threads: config["params"]["salmon"]["threads"] resources: runtime="60m", - ncpus=4, + ncpus=1, mem="32G", container: "docker://combinelab/salmon" @@ -79,7 +86,7 @@ rule export_quantification: ), resources: runtime="30m", - ncpus=2, + ncpus=1, mem="8G", container: "docker://ctglabcnr/tximport" diff --git a/workflow/rules/reads_trimming.smk b/workflow/rules/reads_trimming.smk index 72f84e1..b43b17a 100755 --- a/workflow/rules/reads_trimming.smk +++ b/workflow/rules/reads_trimming.smk @@ -1,26 +1,138 @@ import os -# Althought they may seem equivalent, explicit extra parameters could be defined which will benefit of -# PE or SE sequencing. -rule trimming_pe: +# Determine output filenames based on execution mode +# In CI mode, fastp outputs directly to final names (skip rRNA removal) +# In full mode, fastp outputs to intermediate files for sortmerna processing +if config.get("execution_mode") == "CI": + _trimmed_r1_suffix = "{patient}_1.fastq.gz" + _trimmed_r2_suffix = "{patient}_2.fastq.gz" +else: + _trimmed_r1_suffix = "{patient}_trimmed_1.fastq.gz" + _trimmed_r2_suffix = "{patient}_trimmed_2.fastq.gz" + + +rule trimming: input: - sample=[ - os.path.join(config["resources"]["FASTQ"], "{patient}_1.fastq.gz"), - os.path.join(config["resources"]["FASTQ"], "{patient}_2.fastq.gz"), - ], + unpack(get_fastq), output: - trimmed=[ - os.path.join(config["datadirs"]["trimmed_reads"], "{patient}_1.fastq.gz"), - os.path.join(config["datadirs"]["trimmed_reads"], "{patient}_2.fastq.gz"), - ], - html=os.path.join(config["datadirs"]["trimming_report"], "{patient}_fastp.html"), - json=os.path.join(config["datadirs"]["trimming_report"], "{patient}_fastp.json"), + r1=temp( + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + _trimmed_r1_suffix + )), + r2=temp( + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + _trimmed_r2_suffix + )), + html=os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + "{patient}_fastp.html" + ), + json=os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + "{patient}_fastp.json" + ), params: - extra=config["params"]["fastp"]["pe"], - threads: config["params"]["thread"], + extra=config["params"]["fastp"]["extra"], + threads=config["params"]["fastp"]["threads"], + resources: + mem="20G", + runtime="240m", + ncpus=1, container: "docker://danilotat/eneo" + conda: + "../envs/fastp.yml" log: - os.path.join(config["datadirs"]["logs"]["trimming"], "{patient}.log"), - wrapper: - "v1.0.0/bio/fastp" + os.path.join( + config["datadirs"]["logs"]["trimming"], + "{patient}.log"), + shell: + """ + fastp -i {input.r1} -I {input.r2} \ + -o {output.r1} -O {output.r2} \ + -h {output.html} -j {output.json} \ + -w {params.threads} \ + {params.extra} + """ + + +# rRNA removal step - only included in full mode +if config.get("execution_mode") != "CI": + rule remove_rrna: + input: + r1=os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + "{patient}_trimmed_1.fastq.gz" + ), + r2=os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + "{patient}_trimmed_2.fastq.gz" + ), + rrna_db=config["resources"]["sortmerna_db"], + output: + r1=temp( + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + "{patient}_1.fastq.gz" + )), + r2=temp( + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + "{patient}_2.fastq.gz" + )), + stats=os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + "{patient}_sortmerna.log" + ), + params: + workdir=os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["trimmed_reads"], + "{patient}_sortmerna" + ), + out_prefix=lambda wc, output: os.path.join( + os.path.dirname(output.r1), wc.patient + ), + threads: config["params"]["sortmerna"]["threads"] + resources: + mem="32G", + runtime="240m", + ncpus=1, + container: + "docker://danilotat/sortmerna:latest" + conda: + "../envs/sortmerna.yml" + log: + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["logs"]["trimming"], + "{patient}_sortmerna.log" + ), + shell: + """ + sortmerna \ + --ref {input.rrna_db} \ + --reads {input.r1} \ + --reads {input.r2} \ + --workdir {params.workdir} \ + --aligned {params.workdir}/rrna \ + --other {params.out_prefix} \ + --paired_in \ + --fastx \ + --threads {threads} \ + --out2 2>&1 | tee {output.stats} + mv {params.out_prefix}_fwd.fq.gz {output.r1} + mv {params.out_prefix}_rev.fq.gz {output.r2} + rm -rf {params.workdir} + """ diff --git a/workflow/rules/reporting.smk b/workflow/rules/reporting.smk new file mode 100644 index 0000000..7f1e68e --- /dev/null +++ b/workflow/rules/reporting.smk @@ -0,0 +1,45 @@ +import os +from pathlib import Path + +rule multiqc: + input: + unpack(get_multiqc_inputs), + output: + html=os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["qc_reports"], + "multiqc_report.html" + ), + data=directory( + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["qc_reports"], + "multiqc_data" + ) + ), + params: + outdir=lambda wc, output: Path( + output.data).parent.absolute(), + title="ENEO QC Report", + container: + "docker://ewels/multiqc:latest" + conda: + "../envs/multiqc.yml" + resources: + mem="8G", + runtime="60m", + ncpus=1, + log: + os.path.join( + config["OUTPUT_FOLDER"], + config["datadirs"]["logs"]["trimming"], + "multiqc.log" + ), + shell: + """ + multiqc {input} \ + --outdir {params.outdir} \ + --title "{params.title}" \ + --force \ + 2>&1 | tee {log} + """ diff --git a/workflow/rules/strelka.smk b/workflow/rules/strelka.smk index 49dbe8c..58c8eed 100644 --- a/workflow/rules/strelka.smk +++ b/workflow/rules/strelka.smk @@ -94,7 +94,7 @@ rule Strelka2: ), resources: runtime="240m", - ncpus=2, + ncpus=1, mem="16G", shell: """ @@ -140,7 +140,7 @@ rule SelectStrelka2Calls: ), resources: runtime="20m", - ncpus=2, + ncpus=1, mem="8G", shell: """ diff --git a/workflow/scripts/merge_salmon_quantification.R b/workflow/scripts/merge_salmon_quantification.R index 2ff2c3f..53dcfb3 100755 --- a/workflow/scripts/merge_salmon_quantification.R +++ b/workflow/scripts/merge_salmon_quantification.R @@ -1,39 +1,54 @@ library("dplyr") library("tximport") -library("rtracklayer") +library("rtracklayer") library("tibble") files <- file.path(snakemake@input[["quant"]]) patients <- snakemake@params[["patients"]] gtf_file <- file.path(snakemake@input[["annotation"]]) -coldata <- data.frame(files = files, names = patients, stringsAsFactors = FALSE) + +# Label each quantification file with its patient so the output columns carry +# the sample identity instead of anonymous defaults. +names(files) <- patients gtf <- rtracklayer::import(gtf_file) tx2gene <- as.data.frame(mcols(gtf)[, c("transcript_id", "gene_id")]) - tx2gene <- tx2gene[!is.na(tx2gene$transcript_id), ] -txi <- tximport(files, - type = "salmon", - tx2gene = tx2gene, - ignoreTxVersion = TRUE, - countsFromAbundance = "lengthScaledTPM") -transcript_level_TPM <- txi$abundance %>% +# The GTF carries one row per feature (exon/CDS/...), so each transcript->gene +# mapping is repeated many times; keep the unique pairs only. +tx2gene <- unique(tx2gene) + +# Transcript-level import: txOut = TRUE keeps the per-transcript estimates +# instead of summarising them to gene level. +txi_tx <- tximport(files, + type = "salmon", + txOut = TRUE, + ignoreTxVersion = TRUE, + countsFromAbundance = "lengthScaledTPM") + +transcript_level_TPM <- txi_tx$abundance %>% as.data.frame() %>% tibble::rownames_to_column(var = "transcripts") -write.table(transcript_level_TPM, +write.table(transcript_level_TPM, file = snakemake@output[["transcript"]], - sep = "\t", - row.names = FALSE, + sep = "\t", + row.names = FALSE, quote = FALSE) -gene_level_TPM <- txi$abundance %>% +# Summarise the transcript-level estimates to gene level. +txi_gene <- summarizeToGene(txi_tx, + tx2gene = tx2gene, + ignoreTxVersion = TRUE, + countsFromAbundance = "lengthScaledTPM") + +gene_level_TPM <- txi_gene$abundance %>% as.data.frame() %>% tibble::rownames_to_column(var = "genes") -write.table(gene_level_TPM, +write.table(gene_level_TPM, file = snakemake@output[["gene"]], - sep = "\t", - row.names = FALSE, - quote = FALSE) \ No newline at end of file + sep = "\t", + row.names = FALSE, + quote = FALSE) diff --git a/workflow/supplementary_res/intervals_coding.BED.gz b/workflow/supplementary_res/intervals_coding.BED.gz index 5951a2f..8018983 100644 Binary files a/workflow/supplementary_res/intervals_coding.BED.gz and b/workflow/supplementary_res/intervals_coding.BED.gz differ diff --git a/workflow/supplementary_res/intervals_coding.BED.gz.tbi b/workflow/supplementary_res/intervals_coding.BED.gz.tbi index 9826b60..7b8c61a 100644 Binary files a/workflow/supplementary_res/intervals_coding.BED.gz.tbi and b/workflow/supplementary_res/intervals_coding.BED.gz.tbi differ