Contamination removal, fastq preprocessing

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

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

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

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

.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:

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
+WORKDIR /opt

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)
 

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

pytest.ini

@@ -0,0 +1,3 @@
+[pytest]
+markers =
+    network: marks tests requiring network access (deselect with '-m "not network"')

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

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

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

tests/requirements.txt

@@ -0,0 +1,3 @@
+pytest>=7.0
+pyyaml>=6.0
+requests>=2.28

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)

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