Codex/python conversion

.gitignore

@@ -106,6 +106,22 @@ ENV/
 .idea/
 
 tests/
+!tests/
+tests/*
+!tests/fixtures/
+!tests/fixtures/**
+!tests/test_mini_fixtures_manifest.py
+!tests/test_split_merge_python_expected.py
+!tests/test_parse_se_python.py
+!tests/test_cwl_yaml_adapter.py
+tests/__pycache__/
+tests/*.pyc
+
+.conda-env/
+.snakemake/
+results/
+tests/fixtures/mini/source/failed_jobs_list.txt
+tests/fixtures/mini/source/failed_jobs_list.txt.lck
 
 .ipynb_checkpoints/
 SOURCEME

README.md

@@ -1,111 +1,128 @@
 # repetitive-element-mapping
-pipeline for mapping repetitive elements
-
-# Requirements:
-- Bowtie2=2.2.6
-- Tested with perl=5.22.0 and perl=5.10.1 (perl 5.18+ introduces some randomness in
-iterating over hashes, which may yield slightly different results). This is
-due to the way this pipeline accesses reads via hashes, and determines ties
-in quality by selecting the first read it encounters at the de-duplication
-step. We have tried to mitigate this randomness by iterating over sorted
-hash keys in 0.0.2. See NOTES.
-- cwlref-runner=1.0
-- python 2.7 with pandas and numpy installed
-
-# Installation:
-- For Yeo Lab: ```module load ecliprepmap```
-- For all others:
-    - see the ```DockerRequirement``` from each tool to see the specific requirements for each step. The custom scripts have been tested using the following software/versions:
-      - python=3.6
-      - perl (5.10.1+)
-      - bowtie2=2.2.6
-      - pandas (1.1.3+)
-      - numpy (1.18+)
-      - cwltool==1.0.20180306140409
-      - cwltest
-      - galaxy-lib==17.9.3
-    - ensure that the ```bin/```, ```bin/perl/```, and ```wf/```, and ```cwl/```
-    are correctly in your path.
-
-# Example data: 
-- [example reference data for GRCh38 (hg38)](https://external-collaborator-data.s3-us-west-1.amazonaws.com/reference-data/repeat-family-mapping-grch38.tar.gz)
-- [example input files from eCLIP for GRCh38 (hg38)](https://external-collaborator-data.s3-us-west-1.amazonaws.com/reference-data/example_data_for_repeat_mapping_hg38.tar.gz)
-# Methods:
-- (```map_repetitive_elements_/*.cwl``` - ```parse_bowtie2_output_realtime_includemultifamily_/*.pl```) Runs bowtie2 using the following commands: 
+
+Python + Snakemake pipeline for repetitive element mapping with deterministic mini fixtures for validation.
+
+## What is in this branch
+- Python implementation scripts in `bin/python/`.
+- Snakemake workflow entrypoint in `Snakefile`.
+- Rule modules in `workflow/rules/`.
+- Validation scripts in `workflow/scripts/`.
+- Mini fixtures and expected outputs in `tests/fixtures/mini/`.
+
+## Repository layout
+- `Snakefile`: top-level workflow entrypoint.
+- `config/config.yaml`: runtime file paths for mini validation inputs and expected outputs.
+- `bin/python/split_bam_to_subfiles_SEorPE.py`: split SAM/BAM by UMI prefix.
+- `bin/python/merge_multiple_parsed_files.simplified_20191022.py`: merge parsed statistics files.
+- `workflow/rules/se_foundation.smk`: currently implemented, testable workflow stages.
+- `workflow/scripts/verify_split_manifest.py`: checksum validation for split outputs.
+- `workflow/scripts/verify_merge_against_expected.py`: checksum validation for merged output.
+
+## Requirements
+- Python 3.11+
+- Snakemake 9+
+- samtools 1.23+
+- bowtie2 2.5+
+- pytest 9+
+
+## Installation
+
+### 1) Clone
+```bash
+git clone https://github.com/YeoLab/repetitive-element-mapping.git
+cd repetitive-element-mapping
+git checkout codex/python-conversion
+```
+
+### 2) Create environment
+```bash
+mamba create -y -p ./.conda-env -c conda-forge -c bioconda \
+  python=3.11 snakemake bowtie2 samtools pandas numpy pyyaml pytest
+```
+
+### 3) Activate
+```bash
+mamba activate ./.conda-env
+```
+
+## Input specification
+
+## CWL-style YAML Input Alignment
+This branch supports two equivalent ways to provide run configuration:
+
+1) Populate top-level keys in `config/config.yaml`.
+2) Pass an existing CWL-style job YAML at runtime:
+```bash
+HOME=$(pwd) ./.conda-env/bin/snakemake -j1 -p all \
+  --config cwl_input_yaml=/absolute/path/to/job.yaml
+```
+
+Supported CWL-compatible keys:
+- `barcode1r1FastqGz`, `barcode1rmRepBam`
+- `barcode1Inputr1FastqGz`, `barcode1InputrmRepBam`
+- `bowtie2_db`, `bowtie2_prefix`, `fileListFile1`, `gencodeGTF`, `gencodeTableBrowser`, `repMaskBEDFile`
+- `dataset`, `prefixes`, `se_or_pe`
+
+For `class: File` / `class: Directory` objects, the `path` value is extracted. Relative `path` values are resolved relative to the YAML file location.
+
+The separate `mini_validation` block is only for deterministic fixture tests while conversion is in progress.
+See `docs/config_mapping.md` for the one-to-one mapping details.
+The current implemented Snakemake stages use mini fixture inputs configured in `config/config.yaml`:
+
+- `mini_validation.source_sam_gz`: compressed SAM-like file used for split stage.
+- `mini_validation.merge_input_1`: parsed stats input file #1 for merge stage.
+- `mini_validation.merge_input_2`: parsed stats input file #2 for merge stage.
+- `mini_validation.split_manifest`: expected checksums for 25 split output files.
+- `mini_validation.merge_expected`: expected merged parsed output.
+
+You can update these paths to point to your own test data as long as formats match.
+
+## Output specification
+Running the current workflow stages produces:
+
+- `results/mini/split/python_tmp/*.tmp`: 25 prefix-split files (`AA`..`NN`).
+- `results/mini/split/verified.ok`: split-stage validation success marker.
+- `results/mini/merge/merged.python.parsed`: merged parsed output.
+- `results/mini/merge/verified.ok`: merge-stage validation success marker.
+
+## Run instructions (deployment)
+
+### Local execution
+```bash
+HOME=$(pwd) ./.conda-env/bin/snakemake -j1 -p all
+```
+
+`HOME=$(pwd)` keeps Snakemake cache files inside workspace and avoids host cache permission issues.
+
+### Re-run all stages from scratch
+```bash
+HOME=$(pwd) ./.conda-env/bin/snakemake -j1 -p -F all
+```
+
+### Cluster deployment pattern
+Use your site profile/launcher as usual, for example:
+```bash
+HOME=$(pwd) ./.conda-env/bin/snakemake --profile <your-profile> all
+```
+
+## Testing
+
+### Unit/integration tests
 ```bash
-bowtie2 -q --sensitive -a -p 3 --no-mixed --reorder -x $bowtie_db -1 $fastq_file1 -2 $fastq_file2 2> $bowtie_out
+./.conda-env/bin/python -m pytest -q tests/test_mini_fixtures_manifest.py
 ```
 
-where:
-
-    - $fastq_file1 is read 1 of a trimmed CLIPSEQ expt
-    - $fastq_file2 is read 2 of a trimmed CLIPSEQ expt
-    - $bowtie_db is a bowtie database created using a manually curated set of repeat elements from RepBase or other.
-    - For every proper read pair:
-        - Mismatch score equals the sum of mismatch scores (as determined by the AS: field) for both reads
-        - Reads are compared on the basis of their mismatch scores and read quality, keeping the best or prioritize if equal.
-            - priority is based on input repeat database (e.g. primary transcripts are prioritized over pseudogenes). It is based on the ordering of the repeat database file (MASTER_filelist.wrepbaseandtRNA.enst2id.fixed.UpdatedSimpleRepeat)
-        - Since one read may map to multiple elements, we must compare each element
-            - If read maps to multiple elements of the same family with equal score:
-                - Read mapping information (fastq line) is kept for the element with the highest priority.
-                - Name of all elements is kept
-            - Note: all reads that map to multiple families (not used for downstream analysis).
-        - Create a SAM-like file
-- (```deduplicate.cwl``` - ```duplicate_removal_inline_paired.count_region_other_reads_masksnRNAs_andreparse_SEandPE_20201210_simple.pl```) Merge repeat analysis with unique genomic mapping and remove PCR duplicates
-    - Use the randomer UMIs to remove duplicates based on quality
-        - To save memory, both the SAM-like file and the uniquely mapped BAM file
-        are split based on the first two nucleotides of the randomner umi. De-duplication
-        is performed serially for AA, AC, AG, AT .. NN.
-    - Between a read mapping to both unique genomic and repeat family, if
-    unique mapping to genome is more than 2 mismatches per
-    read = 2 * 2 * 6 alignment score better than to repeat element,
-    throw out repeat element and use genome mapping. Otherwise keep the repeat
-    element mapped read.
-
-### Determining the most correct assignment among elements within one family:
-- Between longer and shorter transcripts: keep the shorter one
-- If a read maps to two places on the same transcript, keep the first
-- Treat "rRNA extra hash" different. See: ```parse_bowtie2_output_realtime_includemultifamily.pl``` $rRNA_extra_hash
-    - This is different because the rRNA precusor transcript contains 18S, 28S, 5.8S transcripts which are treated as separate families.
-# Outputs:
-.parsed file: tabbed file containing 4 comment lines and 4 or 6 columns:
-- #READINFO (AllReads): total number of reads mapped
-- #READINFO (UsableReads): total number of de-duplicated mapped reads
-- #READINFO (GenomicReads): total number of uniquely mapped genomic reads
-- #READINFO (RepFamilyReads): total number of repetitive family mapped reads
-- total_or_element: this is either TOTAL or ELEMENT
-    - Total: family (you can create this by summing all of the elements within this family)
-        - total
-        - family name
-        - number of reads
-        - reads per million
-    - Element: element
-        - element
-        - family name (eg. family1). Can have multiple pipe-delimited families
-        - number of reads
-        - reads per million
-        - family1||transcript1|transcript2|transcript3 (transcript 1/2/3 are all members of family1)
-        - gene name of the transcripts
-.reparsed file: tabbed file with similar structure to .parsed, but without the (Allreads) total.
-.nopipes.tsv: tabbed file containing only non-ambiguously mapped families with the following columns:
-- element: the repeat element
-- IP_read_num: number of reads mapping to the element in IP
-- IP_clip_rpr: number of mapped reads / number of total usable mapped reads
-- Input_read_num: number of reads mapping to the element in Input. +1 pseudocount if no input reads were mapped.
-- Input_clip_rpr: number of mapped input reads / number of total usable mapped input reads
-- Fold_enrichment: IP_clip_rpr / Input_clip_rpr
-- Information_content: IP_clip_rpr * log2(IP_clip_rpr / Input_clip_rpr)
-.withpipes.tsv: tabbed file containing ambiguously and non-ambiguously mapped families
-
-# Notes:
-- Elements (second column) that contain pipes ```|``` are multifamily mapped
-which means that the read is ambiguously mapped to each element. These can generally
-be excluded from downstream analysis.
-- You can sort by information content and log2 fold enrichment.
-Information content is calculated as: log2(clip_rpr/input_rpr).
-- High fold changes and information content (no guidelines yet to this cutoff)
-typically indicate elements enriched in the experiment.
-
-# References:
-- Van Nostrand, E.L., Pratt, G.A., Yee, B.A. et al. Principles of RNA processing from analysis of enhanced CLIP maps for 150 RNA binding proteins. Genome Biol 21, 90 (2020). https://doi.org/10.1186/s13059-020-01982-9
+### Workflow validation run
+```bash
+HOME=$(pwd) ./.conda-env/bin/snakemake -j1 -p -F all
+```
+
+## Regenerating expected artifacts
+If you intentionally update core Python logic, regenerate mini expected artifacts and commit them:
+```bash
+python3 scripts/generate_python_expected_artifacts.py
+```
 
+This refreshes:
+- `tests/fixtures/mini/expected/split.expected.manifest.tsv`
+- `tests/fixtures/mini/expected/merged.expected.parsed`

Snakefile

@@ -0,0 +1,9 @@
+configfile: 'config/config.yaml'
+
+from workflow.config_adapter import merge_cwl_job_yaml_into_config
+
+
+if config.get('cwl_input_yaml'):
+    config = merge_cwl_job_yaml_into_config(config, config['cwl_input_yaml'])
+
+include: 'workflow/rules/se_foundation.smk'