diff --git a/.github/workflows/publish-ghcr.yml b/.github/workflows/publish-ghcr.yml index 8bd7a0d..6470d2a 100644 --- a/.github/workflows/publish-ghcr.yml +++ b/.github/workflows/publish-ghcr.yml @@ -4,6 +4,8 @@ on: push: paths: - "**/Dockerfile" + - "images/druid/complete/**" + - ".github/workflows/publish-ghcr.yml" workflow_dispatch: inputs: dockerfile: @@ -42,7 +44,7 @@ jobs: EVENT_NAME: ${{ github.event_name }} BEFORE_SHA: ${{ github.event.before }} DEFAULT_BRANCH: ${{ github.event.repository.default_branch }} - DOCKERFILE: ${{ inputs.dockerfile }} + REQUESTED_DOCKERFILE: ${{ inputs.dockerfile }} GITHUB_REPOSITORY: ${{ github.repository }} GITHUB_REPOSITORY_OWNER: ${{ github.repository_owner }} GITHUB_SHA: ${{ github.sha }} @@ -50,19 +52,24 @@ jobs: set -euo pipefail if [[ "$EVENT_NAME" == "workflow_dispatch" ]]; then - changed_files="$DOCKERFILE" + changed_files="$REQUESTED_DOCKERFILE" elif [[ "$BEFORE_SHA" =~ ^0+$ ]]; then # A branch's first push has no predecessor SHA. Compare it with # the default branch instead of treating every repository file as # newly added. base_sha="$(git merge-base "origin/$DEFAULT_BRANCH" "$GITHUB_SHA")" - changed_files="$(git diff --name-only --diff-filter=AMR "$base_sha" "$GITHUB_SHA")" + changed_files="$(git diff --name-only --diff-filter=AMRD "$base_sha" "$GITHUB_SHA")" else - changed_files="$(git diff --name-only --diff-filter=AMR "$BEFORE_SHA" "$GITHUB_SHA")" + changed_files="$(git diff --name-only --diff-filter=AMRD "$BEFORE_SHA" "$GITHUB_SHA")" fi mapfile -t dockerfiles < <( - printf '%s\n' "$changed_files" | awk '/(^|\/)Dockerfile$/' + printf '%s\n' "$changed_files" | awk ' + /^images\/druid\/complete\// || $0 == ".github/workflows/publish-ghcr.yml" { + print "images/druid/complete/Dockerfile"; next + } + /(^|\/)Dockerfile$/ { print } + ' | sort -u ) if (( ${#dockerfiles[@]} == 0 )); then @@ -101,6 +108,18 @@ jobs: image_tag="$(printf '%s' "$image_tag" | tr '[:upper:]' '[:lower:]' | sed 's/[^a-z0-9_.-]/-/g')" image="ghcr.io/$owner/$image_name" + if [[ "$context" == "images/druid/complete" ]]; then + # Test the same cached build before publishing either GHCR tag. + docker buildx build --platform linux/amd64 --load \ + --file "$dockerfile" --tag druid:ci "$context" + mkdir -p "$RUNNER_TEMP/druid-test" + docker run --rm --platform linux/amd64 \ + --user "$(id -u):$(id -g)" \ + --volume "$RUNNER_TEMP/druid-test:/work" \ + druid:ci druid smoke-test --work /work \ + 2>&1 | tee "$RUNNER_TEMP/druid-test.log" + fi + echo "Publishing $dockerfile as $image:$image_tag and $image:sha-$revision" docker buildx build \ --file "$dockerfile" \ @@ -114,3 +133,17 @@ jobs: --push \ "$context" done + + - name: Upload DRUID validation logs and results + if: always() + uses: actions/upload-artifact@v4 + with: + name: druid-validation-${{ github.sha }} + if-no-files-found: ignore + retention-days: 14 + path: | + ${{ runner.temp }}/druid-test.log + ${{ runner.temp }}/druid-test/**/Analysis/HEK293_R1/*.csv + ${{ runner.temp }}/druid-test/**/Analysis/HEK293_R1/*.pdf + ${{ runner.temp }}/druid-test/**/Analysis/HEK293_R1/*provenance.json + ${{ runner.temp }}/druid-test/**/Analysis/HEK293_R1/*Info.txt diff --git a/CHANGELOG.md b/CHANGELOG.md index 5eb7ba4..a56d4b6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,9 @@ ### Fixed +- Run DRUID CI validation with the GitHub Actions runner's UID/GID so its + root-owned smoke-test directory no longer prevents artifact upload after a + successful image publication. - Pinned the Bento Tools 2.1.3 base image to RAPIDS Singlecell 0.10.9 so its legacy CUDA 12.6 dependency set is not broken by changes to the `latest` tag. - Updated the Bento Tools 2.1.3 image to cryptography 49.0.0 and its required @@ -17,5 +20,14 @@ ### Added +- Singularity-first DRUID deployment and tutorial instructions for Dockerless + clusters, including the tested `sha-7941ba6aebf0fd4beb48643ec373374b50b02bbc` + image pull command. +- DRUID complete-analysis container at `images/druid/complete`, including + verified GEO/ENA sample mappings, compatibility fixes, and full workflow tests. +- DRUID builds triggered by supporting-file changes; GitHub Actions tests the + container before publishing `ghcr.io/yeolab/druid:complete` and a commit tag, + and retains validation artifacts for 14 days. + - GitHub Actions publishing of changed Dockerfiles to GitHub Container Registry (GHCR), including source, SBOM, and provenance metadata. diff --git a/README.md b/README.md index 9bc5856..d909f81 100644 --- a/README.md +++ b/README.md @@ -70,3 +70,31 @@ apptainer pull fastp_0.23.3.sif docker://ghcr.io/yeolab/fastp:sha-&1 | tee logs/full-analysis.log +``` + +Run on an x86-64 cluster compute node with at least 32 GB RAM (40 GB +recommended) and eight CPUs. The CI test uses a small synthetic reference and +does not download or analyze the full GEO experiment. See the image's +[README](images/druid/complete/README.md) for scheduler guidance and separate +download, reference, analysis, and fitting commands. The image is also tagged +`ghcr.io/yeolab/druid:sha-`. diff --git a/images/druid/complete/.dockerignore b/images/druid/complete/.dockerignore new file mode 100644 index 0000000..800e8d5 --- /dev/null +++ b/images/druid/complete/.dockerignore @@ -0,0 +1,9 @@ +* +!Dockerfile +!scripts/ +!scripts/** +!tests/ +!tests/** +!config/ +!config/** +**/__pycache__ diff --git a/images/druid/complete/.gitignore b/images/druid/complete/.gitignore new file mode 100644 index 0000000..9d100af --- /dev/null +++ b/images/druid/complete/.gitignore @@ -0,0 +1,12 @@ +.DS_Store +upstream/ +reference/ +vendor/ +logs/ +work/ +tmp/ +ExternalData/ +RawSequencingData/ +Analysis/ +__pycache__/ +*.pyc diff --git a/images/druid/complete/CHANGELOG.md b/images/druid/complete/CHANGELOG.md new file mode 100644 index 0000000..910dfce --- /dev/null +++ b/images/druid/complete/CHANGELOG.md @@ -0,0 +1,15 @@ +# Changelog + +## 2026-09-08 + +- Document Singularity-first deployment and tutorial execution on Dockerless + clusters using the tested `sha-7941ba6aebf0fd4beb48643ec373374b50b02bbc` image. + +- Add GitHub Actions publication to `ghcr.io/yeolab/druid:complete` after end-to-end validation, with rebuilds triggered by supporting-file changes. + +- Add an amd64 Dockerfile and Compose service for DRUID's complete-analysis workflow. +- Add verified GEO/ENA mappings and resumable, checksum-verified FASTQ downloads for HEK293 replicate 1. +- Automate combined-reference construction, QC, trimming, alignment, exon/intron quantification, and both half-life analyses. +- Update the GenomicFeatures selection API, port coverage calculation to Python 3, fix exact unique-read filtering, and correct DRUID normalization dimensions while retaining upstream originals. +- Add a synthetic end-to-end test, known-half-life and sequencing-depth-invariance checks, and coverage-statistics checks. +- Document resource requirements, compatibility changes, provenance, and limits of validation. diff --git a/images/druid/complete/Dockerfile b/images/druid/complete/Dockerfile new file mode 100644 index 0000000..1602b2c --- /dev/null +++ b/images/druid/complete/Dockerfile @@ -0,0 +1,47 @@ +# syntax=docker/dockerfile:1 +# UCSC's distributed tools are x86-64. Build with --platform linux/amd64. +FROM debian:bookworm-slim@sha256:88200866dfff7ea7f5cbcb6ec7c8a701889efe6fe859fe64d6990e4b07ea4171 +ARG TARGETARCH +ARG DRUID_COMMIT=455be347cab87b62a2d1f081ec34cafffdf81622 +ENV DEBIAN_FRONTEND=noninteractive LC_ALL=C.UTF-8 LANG=C.UTF-8 \ + DRUID_HOME=/opt/DRUID PATH=/opt/DRUID:/opt/DRUID/Scripts:/opt/workflow:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin \ + THREADS=8 TMPDIR=/work/tmp +SHELL ["/bin/bash", "-o", "pipefail", "-c"] +RUN test "$TARGETARCH" = amd64 || { echo 'Build with --platform linux/amd64'; exit 1; } +RUN mkdir -p /work/tmp && apt-get update && apt-get install -y --no-install-recommends \ + ca-certificates curl wget git perl bash coreutils gawk parallel procps \ + samtools bedtools rna-star fastqc trimmomatic default-jre-headless \ + python3 python3-numpy python3-htseq python-is-python3 sra-toolkit \ + r-base-core r-cran-gplots r-cran-plyr r-bioc-genomicfeatures r-bioc-rtracklayer \ + libmariadb3 libpng16-16 libssl3 libuuid1 libgomp1 \ + && dpkg-query -W > /opt/debian-packages.tsv \ + && rm -rf /var/lib/apt/lists/* +RUN git init /opt/DRUID && cd /opt/DRUID \ + && git remote add origin https://github.com/risslandlab/DRUID.git \ + && git fetch --depth 1 origin "$DRUID_COMMIT" && git checkout --detach FETCH_HEAD \ + && test "$(git rev-parse HEAD)" = "$DRUID_COMMIT" \ + && echo "$DRUID_COMMIT" > /opt/DRUID/UPSTREAM_COMMIT && rm -rf .git +# Checksums intentionally fail if UCSC replaces these unversioned downloads. +RUN curl -fL --retry 5 https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa -o /usr/local/bin/twoBitToFa \ + && curl -fL --retry 5 https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/genePredToGtf -o /usr/local/bin/genePredToGtf \ + && echo '10b2cd0d1f8208e8c0e9990b4a97d90cda4eec754b13930a10d55a40bda56710 /usr/local/bin/twoBitToFa' | sha256sum -c - \ + && echo '3596b456da9c423686a2a3f01edab9915d1969144601c0aac3c8b694375a2c18 /usr/local/bin/genePredToGtf' | sha256sum -c - \ + && chmod 755 /usr/local/bin/twoBitToFa /usr/local/bin/genePredToGtf +COPY scripts/ /opt/workflow/ +COPY tests/ /opt/tests/ +COPY config/ /opt/config/ +RUN python3 /opt/workflow/patch_upstream.py \ + && chmod +x /opt/DRUID/*.sh /opt/DRUID/*.py /opt/DRUID/Scripts/* /opt/workflow/*.sh \ + && ln -s /usr/share/java/trimmomatic.jar /opt/trimmomatic.jar \ + && ln -s /usr/share/trimmomatic /opt/trimmomatic \ + && ln -s /opt/workflow/druid.py /usr/local/bin/druid \ + && chmod +x /opt/workflow/druid.py \ + && Rscript -e 'library(GenomicFeatures); library(rtracklayer); library(plyr); library(gplots); source("/opt/DRUID/DRUID.R"); stopifnot(is.function(DRUID))' \ + && python3 -c 'import numpy, HTSeq' \ + && STAR --version && samtools --version && bedtools --version && fastqc --version \ + && mkdir -p /work/tmp +LABEL org.opencontainers.image.title="DRUID complete analysis" \ + org.opencontainers.image.source="https://github.com/risslandlab/DRUID" \ + org.opencontainers.image.revision="455be347cab87b62a2d1f081ec34cafffdf81622" +WORKDIR /work +CMD ["bash"] diff --git a/images/druid/complete/LICENSE-DRUID b/images/druid/complete/LICENSE-DRUID new file mode 100644 index 0000000..a4471a5 --- /dev/null +++ b/images/druid/complete/LICENSE-DRUID @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2017 risslandlab + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/images/druid/complete/README.md b/images/druid/complete/README.md new file mode 100644 index 0000000..7e86e04 --- /dev/null +++ b/images/druid/complete/README.md @@ -0,0 +1,171 @@ +# DRUID complete-analysis image + +Singularity-ready command-line image for the [DRUID complete-analysis tutorial](https://github.com/risslandlab/DRUID/wiki/Complete-analysis), pinned to upstream commit `455be347cab87b62a2d1f081ec34cafffdf81622`. + +This guide assumes a cluster where Docker is unavailable. Singularity runs the +published OCI image without a Docker daemon. Your project directory is bound +at `/work`; input data, references, temporary files, and analysis outputs stay +in that directory. The image contains software only. + +## Published image + +GitHub Actions builds this directory from `YeoLab/containers` and publishes +`ghcr.io/yeolab/druid:complete` plus `ghcr.io/yeolab/druid:sha-` +after the end-to-end tests pass. All supporting-file changes trigger CI. Pull +the tested revision used by this guide on a login node or any node with network +access: + +```bash +singularity pull druid.sif docker://ghcr.io/yeolab/druid:sha-7941ba6aebf0fd4beb48643ec373374b50b02bbc +export DRUID_IMAGE="$PWD/druid.sif" +``` + +`druid.sif` is an immutable local image file. Keep it outside the analysis +directory if your cluster backs that directory up or applies a quota to it. To +build the image from this repository instead, use `docker compose build` inside +`images/druid/complete/` on a Docker-capable development machine. + +If GHCR returns an authorization error, authenticate Singularity with a GitHub +token that has `read:packages` access, or ask a package administrator to make +the package public. The repository [GHCR authentication instructions](../../../README.md#pulling-a-singularityapptainer-image-from-ghcr) +cover the login command. + +## Test the pulled image + +```bash +mkdir -p "$HOME/druid-smoke-test" +# Includes a small synthetic reference and all six time points, through half-life fitting. +singularity exec --cleanenv \ + --bind "$HOME/druid-smoke-test:/work" --pwd /work \ + "$DRUID_IMAGE" druid smoke-test --work /work +``` + +The image targets **linux/amd64**, because the bundled UCSC binaries are +x86-64. Run it on an x86-64 Linux cluster node. The smoke-test outputs are +written under `$HOME/druid-smoke-test` and can be removed after inspection. +For maintainers with Docker, a direct build is also supported: + +```bash +docker build --platform linux/amd64 -t druid:complete . +``` + +## Run the tutorial data + +Request **at least 32 GB RAM; 40 GB is recommended**, with eight CPUs from +your scheduler. The full human/fly/yeast workflow checks the effective memory +limit before running. Plan for roughly **100 GB of free working space**, +depending on BAM sizes and temporary coverage files. The six compressed FASTQs +alone total 1,615,254,957 bytes (~1.50 GiB). + +Run on an allocated compute node, not a shared login node. For Slurm, a typical +interactive allocation is: + +```bash +srun --pty --cpus-per-task=8 --mem=40G --time=24:00:00 bash +``` + +Then create or enter the directory that will hold the tutorial inputs and +outputs. Set `THREADS` to no more than the CPUs requested from the scheduler. + +```bash +mkdir -p "$HOME/DRUID" +cd "$HOME/DRUID" +mkdir -p tmp logs +# Full workflow: FASTQ download -> reference construction -> QC -> trimming -> +# alignment -> HTSeq -> preDRUID -> DRUID and spike-in half-life calculations. +set -o pipefail +THREADS=8 singularity exec --cleanenv \ + --env THREADS=8,TMPDIR=/work/tmp \ + --bind "$PWD:/work" --pwd /work \ + "$DRUID_IMAGE" druid all 2>&1 | tee logs/full-analysis.log +``` + +The same work can be run in separate stages: + +```bash +# Can run before requesting the full-memory compute allocation. +singularity exec --cleanenv --env THREADS=8,TMPDIR=/work/tmp \ + --bind "$PWD:/work" --pwd /work "$DRUID_IMAGE" druid download +# Require the full allocation for the remaining stages. +singularity exec --cleanenv --env THREADS=8,TMPDIR=/work/tmp \ + --bind "$PWD:/work" --pwd /work "$DRUID_IMAGE" druid reference +singularity exec --cleanenv --env THREADS=8,TMPDIR=/work/tmp \ + --bind "$PWD:/work" --pwd /work "$DRUID_IMAGE" druid analyze +# Repeat only the final R analysis from existing count and coverage tables. +singularity exec --cleanenv --env THREADS=8,TMPDIR=/work/tmp \ + --bind "$PWD:/work" --pwd /work "$DRUID_IMAGE" druid fit +``` + +`DRUID_IMAGE` records the absolute path to the file created by `singularity +pull`; use the same bind and working-directory arguments for every DRUID +command. +`--cleanenv` avoids accidentally inheriting host Python, R, or Java settings. +Samples are processed sequentially to bound memory; FastQC uses at most four +simultaneous jobs. All data are single-end, reverse-stranded, as in the +tutorial. QC reports still need scientific review. + +The downloader verifies ENA MD5 checksums, skips verified existing files, and resumes `.partial` downloads. A checksum failure stops processing. If a partial file has a checksum mismatch, remove that partial file and retry. Analysis refuses to overwrite an existing `Analysis/HEK293_R1` directory; after a failed run, inspect its logs and move the directory aside before rerunning `analyze`. Reference preparation may be rerun, using cached downloads. `all` rebuilds the reference and expects a fresh analysis directory. + +## Samples and reference + +`config/HEK293_R1.tsv` maps the six GEO samples to ENA-hosted FASTQs from the same SRA study, **SRP108376 / PRJNA388686**. The GEO family metadata was used to verify HEK293 cells, 4SU selection, replicate 1, and treatment times. Steady-state, replicate 2, inhibitor-treated, mouse, and abnormal-spike-in samples are excluded from this tutorial run. + +| Hours | Local sample | GEO | SRA run | +|---:|---|---|---| +| 1 | HEK293_R1_01 | GSM2645366 | SRR5634082 | +| 2 | HEK293_R1_02 | GSM2645367 | SRR5634083 | +| 4 | HEK293_R1_04 | GSM2645368 | SRR5634084 | +| 8 | HEK293_R1_08 | GSM2645369 | SRR5634085 | +| 12 | HEK293_R1_12 | GSM2645370 | SRR5634086 | +| 24 | HEK293_R1_24 | GSM2645371 | SRR5634087 | + +Reference preparation downloads UCSC **hg38/refGene, dm6/refGene, and sacCer3/sgdGene**, filters canonical human/fly chromosomes, labels chromosome and feature identifiers with `h`, `f`, or `y`, merges the references, runs `process_GTF.R`, and builds a STAR index with `sjdbOverhang=50` and `genomeSAsparseD=2`. UCSC genePred tables are downloaded over HTTPS and passed to `genePredToGtf` as local files, avoiding the tutorial's remote MySQL dependency. The `bin` column is removed; the yeast-only trailing `proteinID` field is also removed to produce valid genePred input. All downloaded inputs remain under `ExternalData/Genomes`. + +The final FASTA and GTF SHA-256 hashes are recorded in `provenance.json`. UCSC annotations at these URLs can change. To use archived, already filtered and species-labeled references: + +```bash +singularity exec --cleanenv --env THREADS=8,TMPDIR=/work/tmp \ + --bind "$PWD:/work" --pwd /work "$DRUID_IMAGE" druid reference \ + --fasta /work/my-reference/combined.fa --gtf /work/my-reference/combined.gtf +``` + +Do not use `--small-genome` for human references; it changes STAR index parameters specifically for tiny test genomes. The custom reference must retain the tutorial's chromosome/species suffix and exon identifier conventions. + +## Outputs + +`Analysis/HEK293_R1/` contains: + +- `FastQC/`, `TrimmedReads/`, `STAR/`, indexed `BAM/`, and `Logs/Logs.tsv`. +- `htseq_exons_introns/combined_exon.txt` and `combined_intron.txt`. +- `intersect/`: four coverage CSVs per sample. +- `half-lives.DRUID.csv` and `half-lives.spikeins.csv`: standard and leave-one-time-point-out estimates, in hours. Failed or negative fits can be `NA`. +- `DRUID.replicate.1.heatmap.DRUID.pdf` and `.filterIntrons.norm.csv`. +- `analysis.rds`: ordered samples, counts, coverage tables, and results. +- R session information, Debian package versions, sample manifest, and reference hashes. + +The runner explicitly orders all matrices and STAR library sizes by numeric sampling time. It uses the wiki's trimming, alignment, counting, intron-coverage threshold, four-cluster selection, and k-means seed 42. As in upstream, DRUID fitting uses no cell-cycle correction; the spike-in comparison uses a 24-hour doubling time. + +## Compatibility changes + +The container retains original core files as `*.upstream.*` under `/opt/DRUID` and applies these documented changes: + +1. Port `coverage.py` from Python 2 to Python 3, retaining per-base means, population standard deviations, feature lengths, and covered fractions. Accept quoted GTF attributes. +2. Replace the `preDRUID.sh` wrapper's relative paths and shared `/dev/shm` filenames with installed paths and isolated temporary directories under `/work/tmp`. Errors propagate, and unique mapping matches the exact `NH:i:1` tag, excluding `NH:i:10` and higher. +3. Remove the non-executable example section from the sourced `DRUID.R`. +4. Correct DRUID's normalization dimensions: divide each sample column by its library size and each intron row by its maximum. Upstream constructs an intron-by-intron matrix instead of an intron-by-time-point matrix and uses vector recycling for library-size division. This is a numerical correction, so results should not be represented as byte-identical outputs of unpatched upstream. +5. Update `process_GTF.R` to use the current GenomicFeatures `filter` argument in place of the removed `vals` argument, preserving longest-transcript selection; replace the removed `ignoreSelf` overlap argument by explicit diagonal-hit removal. +6. Use a command-line driver in place of manual wiki steps and hard-coded batch helpers. STAR uses `NoSharedMemory`, a 4 GB BAM-sort buffer, and disk temporary files for container compatibility. The biological filtering parameters remain as documented. + +The image uses Debian Bookworm's R/Bioconductor, STAR, HTSeq, samtools, bedtools, FastQC, Trimmomatic, GNU Parallel, and SRA Toolkit. It pins the base image digest, DRUID commit, and UCSC binary checksums. Debian package updates remain possible on rebuild; `/opt/debian-packages.tsv` records exact installed versions. A UCSC download changing will fail its checksum verification rather than silently replacing the binary. + +This is a runnable implementation of the tutorial with modern dependencies, not a claim of exact reproduction of the 2018 paper. The paper used a modified STAR 2.5.2a and historical annotations; its methods also describe manual intron-cluster selection, while the public script/tutorial selects a cluster automatically. The full published-data run requires the memory allocation above. See `VALIDATION.md` for the actual checks performed. + +## Sources + +- [DRUID source and MIT license](https://github.com/risslandlab/DRUID) +- [Complete-analysis tutorial](https://github.com/risslandlab/DRUID/wiki/Complete-analysis) +- [GEO GSE99517](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE99517) +- [ENA study SRP108376](https://www.ebi.ac.uk/ena/browser/view/SRP108376) +- [Lugowski et al., DRUID paper](https://pmc.ncbi.nlm.nih.gov/articles/PMC5900561/) + +Upstream-derived code is covered by `LICENSE-DRUID`. Other included tools retain their respective licenses. diff --git a/images/druid/complete/VALIDATION.md b/images/druid/complete/VALIDATION.md new file mode 100644 index 0000000..f8dda4c --- /dev/null +++ b/images/druid/complete/VALIDATION.md @@ -0,0 +1,57 @@ +# Validation — 2026-09-08 + +The image was built and tested locally using Docker Desktop on Apple Silicon with `linux/amd64` emulation, eight available CPUs, and approximately 12 GB allocated to Docker. + +## Passed + +- `docker compose build`: built `druid:complete` successfully. +- `docker buildx build --check --platform linux/amd64 .`: no warnings. +- `docker compose config --quiet`: valid Compose configuration. +- Python compilation and Bash syntax checks. +- All required R libraries and Python modules load inside the image. +- Coverage statistics match a hand-calculated four-base example, including zero-covered bases and BED/GTF coordinate conversion. +- A BAM fixture verifies exact `NH:i:1` filtering, exclusion of `NH:i:10`, and reverse-stranded coverage. +- Annotation fixtures verify longest-transcript selection, removal of overlapping exons, exclusion of introns overlapping annotated exons, and preservation of opposite-strand features. +- The six-time-point end-to-end test runs GTF processing, STAR indexing, FastQC, Trimmomatic, STAR alignment, samtools validation/indexing, exon/intron HTSeq counting, count merging, preDRUID, DRUID fitting, spike-in fitting, and heatmap generation. +- End-to-end results contain 48 human genes, seven fit columns (standard plus six leave-one-out fits), and 60 introns across the three species. All 48 standard DRUID and spike-in fits are finite. +- Fitting an analytical bounded-growth curve recovers a known three-hour half-life within `1e-4` hours. +- Scaling the six libraries by different factors preserves DRUID half-lives within tolerance `1e-5`. +- The complete current hg38/refGene, dm6/refGene, and sacCer3/sgdGene annotations were merged and processed successfully: 306,488 representative exons, 294,691 non-overlapping exons, and 236,662 filtered introns. This tests the real annotations independently of the unrun full-genome index. See `logs/real-annotations.log` and `reference/validate-real-annotations.py`. +- Both UCSC converters execute on real sacCer3 inputs. The resulting labeled yeast annotation passes `process_GTF.R` and produces nonempty intron output. +- The full-reference command stops with an explicit memory error under the current Docker allocation, before downloading a genome or launching the human index build. + +Successful end-to-end outputs are in `work/smoke.1FV5nN/Analysis/HEK293_R1/`. Logs are in `logs/build.log`, `logs/smoke-test.log`, `logs/dockerfile-check.log`, and `logs/memory-guard.log`. Earlier failed test directories remain under `work/` for inspection and are not production analyses. + +Final tested image ID: + +```text +sha256:93c16c51313f0b81fe05118658ea447168835a9bcb88c1d6ccb4dfabd7ba684f +``` + +## Key installed versions + +| Component | Version | +|---|---| +| DRUID upstream | `455be347cab87b62a2d1f081ec34cafffdf81622`, with documented patches | +| STAR | 2.7.10b | +| samtools | 1.16.1 | +| bedtools | 2.30.0 | +| FastQC | 0.11.9 | +| Trimmomatic | 0.39 | +| Python | 3.11.2 | +| NumPy | 1.24.2 | +| HTSeq | 1.99.2 | +| R | 4.2.2 | +| GenomicFeatures | 1.50.4 | +| rtracklayer | 1.58.0 | +| gplots | 3.1.3 | +| plyr | 1.8.8 | +| SRA Toolkit | 3.0.3 | + +The complete package inventory is `/opt/debian-packages.tsv` in the image and `software-versions.tsv` in each completed analysis. Per-analysis workflow provenance includes the upstream commit and hashes of the executed scripts. + +## Limits + +A complete GSE99517 FASTQ-to-half-life run has **not** been performed. Docker currently has about 12 GB RAM; full indexing and analysis require at least 32 GB. Full compressed FASTQ downloads were also not completed: the slow ENA transfer was stopped, leaving a resumable `RawSequencingData/HEK293_R1/HEK293_R1_01.fastq.gz.partial` file. `druid download` resumes it and verifies its MD5 before using it. + +Synthetic execution and numerical checks establish that the workflow connects and executes, but they do not establish agreement with the paper's published half-lives. Modern package versions, the normalization correction, reference annotation changes, and automatic cluster selection can change numerical results. No human genome mapping, biological reproduction, or full-scale performance claim is made. diff --git a/images/druid/complete/compose.yaml b/images/druid/complete/compose.yaml new file mode 100644 index 0000000..0ba866b --- /dev/null +++ b/images/druid/complete/compose.yaml @@ -0,0 +1,14 @@ +services: + druid: + image: druid:complete + platform: linux/amd64 + build: + context: . + working_dir: /work + volumes: + - .:/work + environment: + THREADS: ${THREADS:-8} + TMPDIR: /work/tmp + init: true + command: ["bash"] diff --git a/images/druid/complete/config/HEK293_R1.tsv b/images/druid/complete/config/HEK293_R1.tsv new file mode 100644 index 0000000..1f514f6 --- /dev/null +++ b/images/druid/complete/config/HEK293_R1.tsv @@ -0,0 +1,7 @@ +sample time_hours geo_sample run url md5 bytes +HEK293_R1_01 1 GSM2645366 SRR5634082 https://ftp.sra.ebi.ac.uk/vol1/fastq/SRR563/002/SRR5634082/SRR5634082.fastq.gz 4b3dfe9ca0e604b458e52e62ab851e6a 308022137 +HEK293_R1_02 2 GSM2645367 SRR5634083 https://ftp.sra.ebi.ac.uk/vol1/fastq/SRR563/003/SRR5634083/SRR5634083.fastq.gz d45cfd3891768dbde534c5ba3a6d26ad 235831929 +HEK293_R1_04 4 GSM2645368 SRR5634084 https://ftp.sra.ebi.ac.uk/vol1/fastq/SRR563/004/SRR5634084/SRR5634084.fastq.gz b3d10d2f4acaa22c8aae467d6ae98bf7 247066051 +HEK293_R1_08 8 GSM2645369 SRR5634085 https://ftp.sra.ebi.ac.uk/vol1/fastq/SRR563/005/SRR5634085/SRR5634085.fastq.gz 4d21f7bbf1defd7a4e7f957c14bc1a29 192884586 +HEK293_R1_12 12 GSM2645370 SRR5634086 https://ftp.sra.ebi.ac.uk/vol1/fastq/SRR563/006/SRR5634086/SRR5634086.fastq.gz fa736ef7066556a65b96a36cecc3d054 303331201 +HEK293_R1_24 24 GSM2645371 SRR5634087 https://ftp.sra.ebi.ac.uk/vol1/fastq/SRR563/007/SRR5634087/SRR5634087.fastq.gz e8d863101e2dad79d979fe72662d4207 328119053 diff --git a/images/druid/complete/scripts/coverage.py b/images/druid/complete/scripts/coverage.py new file mode 100644 index 0000000..d24c2b1 --- /dev/null +++ b/images/druid/complete/scripts/coverage.py @@ -0,0 +1,32 @@ +#!/usr/bin/env python3 +"""Python 3 port of DRUID coverage.py, accepting quoted GTF attributes.""" +import csv +import re +import sys +import numpy as np + +features, info = {}, {} +with open(sys.argv[1], newline='') as handle: + for row in csv.reader(handle, delimiter='\t'): + chrom = re.fullmatch(r'(chr[0-9A-Z]+)([hfy])', row[0]) + if not chrom: + raise ValueError(f'Expected species-suffixed canonical chromosome: {row[0]}') + attrs = dict(re.findall(r'(\w+)\s+"?([^";]+)"?;', row[8])) + gene = attrs['gene_id'][:-1] + exon_id = re.fullmatch(r'(.+)\.(\d+)[hfy]', attrs['exon_id']) + if not exon_id: + raise ValueError(f'Expected transcript.number exon_id: {row[8]}') + transcript, exon = exon_id.groups() + key = (gene, chrom[2], exon) + start, end = int(row[3]) - 1, int(row[4]) + if key not in features: + features[key] = np.zeros(end - start, dtype=np.int32) + info[key] = [chrom[1], chrom[2], gene, transcript, exon] + a, b = max(start, int(row[10])) - start, min(end, int(row[11])) - start + if b - a != int(row[13]): + raise ValueError(f'Intersection length mismatch: {row}') + features[key][a:b] += int(row[12]) +out = csv.writer(sys.stdout, lineterminator='\n') +out.writerow(['chromosome', 'organism', 'gene', 'transcript', 'exon', 'mean', 'stdev', 'size', 'coverage']) +for key, values in features.items(): + out.writerow(info[key] + [np.mean(values), np.std(values), values.size, np.count_nonzero(values) / values.size]) diff --git a/images/druid/complete/scripts/druid.py b/images/druid/complete/scripts/druid.py new file mode 100644 index 0000000..c4260b7 --- /dev/null +++ b/images/druid/complete/scripts/druid.py @@ -0,0 +1,281 @@ +#!/usr/bin/env python3 +"""Reproducible command-line runner for the DRUID complete-analysis tutorial.""" +import argparse +import csv +import hashlib +import json +import os +from pathlib import Path +import re +import shutil +import subprocess +import sys +from datetime import datetime, timezone + +HOME = Path(os.environ.get('DRUID_HOME', '/opt/DRUID')) +HERE = Path(__file__).resolve().parent + + +def run(args, cwd=None, output=None, env=None): + args = [str(a) for a in args] + print('+ ' + ' '.join(args), flush=True) + if output: + with Path(output).open('w') as handle: + subprocess.run(args, cwd=cwd, env=env, stdout=handle, check=True) + else: + subprocess.run(args, cwd=cwd, env=env, check=True) + + +def digest(path, kind='sha256'): + h = hashlib.new(kind) + with Path(path).open('rb') as handle: + for chunk in iter(lambda: handle.read(1024 * 1024), b''): + h.update(chunk) + return h.hexdigest() + + +def download(url, path, md5=None): + path = Path(path) + path.parent.mkdir(parents=True, exist_ok=True) + if path.exists(): + if md5 and digest(path, 'md5') != md5: + raise ValueError(f'Checksum mismatch for existing file: {path}; move it aside before retrying') + print(f'Using {path}', flush=True) + return + partial = path.with_name(path.name + '.partial') + run(['curl', '-fL', '--retry', '5', '--retry-delay', '3', '-C', '-', url, '-o', partial]) + if md5 and digest(partial, 'md5') != md5: + raise ValueError(f'Checksum mismatch: {partial}; remove this partial download before retrying') + partial.rename(path) + + +def samples(args): + with args.manifest.open() as handle: + rows = list(csv.DictReader(handle, delimiter='\t')) + if len(rows) < 4: + raise ValueError('At least four time points are required for leave-one-out fitting') + for r in rows: + if not re.fullmatch(r'[A-Za-z0-9_.-]+', r['sample']): + raise ValueError(f'Invalid sample identifier: {r["sample"]}') + if float(r['time_hours']) <= 0: + raise ValueError('Approach-to-equilibrium time points must be positive') + if len({r['sample'] for r in rows}) != len(rows) or len({float(r['time_hours']) for r in rows}) != len(rows): + raise ValueError('Samples and time points must be unique') + return sorted(rows, key=lambda r: float(r['time_hours'])) + + +def memory_check(): + # Check both Linux physical memory and container cgroup limits. + limits = [] + meminfo = Path('/proc/meminfo') + if meminfo.exists(): + limits.append(int(re.search(r'MemTotal:\s+(\d+)', meminfo.read_text())[1]) * 1024) + for p in ('/sys/fs/cgroup/memory.max', '/sys/fs/cgroup/memory/memory.limit_in_bytes'): + if Path(p).exists(): + value = Path(p).read_text().strip() + if value.isdigit(): + limits.append(int(value)) + if limits and min(limits) < 30 * 1024**3: + raise RuntimeError('The combined human/fly/yeast workflow requires at least 32 GB assigned to Docker (40 GB recommended). The current memory limit is too low.') + + +def fetch_fastq(args): + for r in samples(args): + download(r['url'], args.work / 'RawSequencingData' / args.dataset / (r['sample'] + '.fastq.gz'), r['md5']) + + +def unpack_genepred(annotation, target, table): + import gzip + with gzip.open(annotation, 'rt') as src, open(target, 'w') as dst: + for line in src: + fields = line.rstrip('\n').split('\t')[1:] # discard UCSC bin + # sgdGene appends proteinID, not genePredExt's score/name2 fields. + if table == 'sgdGene': + fields = fields[:10] + dst.write('\t'.join(fields) + '\n') + + +def label_reference(fasta, gtf, suffix, filtered, outfa, outgtf): + # Exactly the tutorial's canonical-chromosome filter; yeast retains all contigs. + keep = False + with open(fasta) as src, open(outfa, 'w') as dst: + for line in src: + if line.startswith('>'): + name = line[1:].strip() + keep = not filtered or bool(re.fullmatch(r'chr[0-9A-Z]{1,2}', name)) + if keep: + dst.write('>' + name + suffix + '\n') + elif keep: + dst.write(line) + with open(gtf) as src, open(outgtf, 'w') as dst: + for line in src: + if line.startswith('#'): + continue + fields = line.rstrip('\n').split('\t') + if len(fields) != 9: + continue + if filtered and not re.fullmatch(r'chr[0-9A-Z]{1,2}', fields[0]): + continue + fields[0] += suffix + fields[8] = re.sub(r'((?:gene_id|transcript_id|exon_id|gene_name) ")([^"]+)(")', + lambda m: m[1] + m[2] + suffix + m[3], fields[8]) + dst.write('\t'.join(fields) + '\n') + + +def prepare_reference(args): + if not args.small_genome: + memory_check() + base = args.work / 'ExternalData' + ref = base / 'Genomes' / 'hg38dm6sacCer3' + index = base / 'IndexedGenomes' / 'hg38dm6sacCer3_50' + ref.mkdir(parents=True, exist_ok=True) + index.mkdir(parents=True, exist_ok=True) + fa, gtf = ref / 'hg38dm6sacCer3.fa', ref / 'hg38dm6sacCer3.gtf' + if bool(args.fasta) != bool(args.gtf): + raise ValueError('Supply both --fasta and --gtf (already labeled and filtered)') + if args.fasta: + shutil.copyfile(args.fasta, fa) + shutil.copyfile(args.gtf, gtf) + else: + fastas, gtfs = [], [] + for genome, table, suffix, filtered in [('hg38', 'refGene', 'h', True), ('dm6', 'refGene', 'f', True), ('sacCer3', 'sgdGene', 'y', False)]: + folder = base / 'Genomes' / genome + folder.mkdir(parents=True, exist_ok=True) + bit = folder / f'{genome}.2bit' + annotation = folder / f'{table}.txt.gz' + download(f'https://hgdownload.soe.ucsc.edu/goldenPath/{genome}/bigZips/{genome}.2bit', bit) + download(f'https://hgdownload.soe.ucsc.edu/goldenPath/{genome}/database/{table}.txt.gz', annotation) + run(['twoBitToFa', bit, folder / 'genome.fa']) + # UCSC table dumps contain a leading bin field; file input avoids remote MySQL. + unpack_genepred(annotation, folder / 'genes.gp', table) + run(['genePredToGtf', '-utr', 'file', folder / 'genes.gp', folder / 'genes.gtf']) + label_reference(folder / 'genome.fa', folder / 'genes.gtf', suffix, filtered, + folder / 'labeled.fa', folder / 'labeled.gtf') + fastas.append(folder / 'labeled.fa') + gtfs.append(folder / 'labeled.gtf') + for sources, target in [(fastas, fa), (gtfs, gtf)]: + with target.open('wb') as out: + for source in sources: + with source.open('rb') as src: + shutil.copyfileobj(src, out) + run(['Rscript', HOME / 'process_GTF.R', gtf]) + for suffix in ('exons.gtf', 'exons.no.gtf', 'introns.gtf'): + if not (ref / f'hg38dm6sacCer3.{suffix}').stat().st_size: + raise ValueError(f'Empty annotation: {suffix}') + run(['samtools', 'faidx', fa]) + with open(str(fa) + '.fai') as handle: + sizes = sorted((line.split('\t')[:2] for line in handle), key=lambda r: r[0]) + (ref / 'sizes.genome').write_text(''.join('\t'.join(r) + '\n' for r in sizes)) + command = ['STAR', '--runThreadN', args.threads, '--runMode', 'genomeGenerate', '--genomeDir', index, + '--genomeFastaFiles', fa, '--sjdbGTFfile', gtf, '--genomeSAsparseD', '2', '--sjdbOverhang', '50', + '--limitGenomeGenerateRAM', '31000000000'] + if args.small_genome: + command += ['--genomeSAindexNbases', '6', '--genomeChrBinNbits', '12'] + run(command, cwd=index) + provenance = {'created_utc': datetime.now(timezone.utc).isoformat(), 'small_genome': args.small_genome, + 'fasta_sha256': digest(fa), 'gtf_sha256': digest(gtf)} + (ref / 'provenance.json').write_text(json.dumps(provenance, indent=2) + '\n') + + +def analyze(args): + rows = samples(args) + ref = args.work / 'ExternalData/Genomes/hg38dm6sacCer3' + index = args.work / 'ExternalData/IndexedGenomes/hg38dm6sacCer3_50' + provenance = json.loads((ref / 'provenance.json').read_text()) + if not provenance['small_genome']: + memory_check() + analysis = args.work / 'Analysis' / args.dataset + if analysis.exists(): + raise RuntimeError(f'Analysis directory already exists: {analysis}. Use `druid fit` to repeat fitting, or move the analysis aside before rerunning.') + raw = args.work / 'RawSequencingData' / args.dataset + fastqs = [raw / (r['sample'] + '.fastq.gz') for r in rows] + for path in fastqs: + if not path.is_file(): + raise FileNotFoundError(path) + for name in ('FastQC', 'TrimmedReads', 'STAR', 'BAM', 'Logs', 'htseq_exons_introns'): + (analysis / name).mkdir(parents=True, exist_ok=True) + (analysis / 'data.files.txt').write_text(''.join(str(p) + '\n' for p in fastqs)) + # Limit simultaneous QC memory; other tools use the full thread budget. + run(['fastqc', '--outdir', analysis / 'FastQC', '--threads', min(args.threads, 4), *fastqs]) + logs = [] + for r, fq in zip(rows, fastqs): + sample = r['sample'] + trimmed = analysis / 'TrimmedReads' / fq.name + run(['java', '-jar', '/opt/trimmomatic.jar', 'SE', '-threads', args.threads, '-phred33', fq, trimmed, + 'ILLUMINACLIP:/opt/trimmomatic/TruSeq3-SE.fa:2:30:10', 'LEADING:3', 'TRAILING:3', 'SLIDINGWINDOW:4:15', 'MINLEN:36']) + star = analysis / 'STAR' / sample + star.mkdir() + run(['STAR', '--runThreadN', args.threads, '--genomeDir', index, '--readFilesIn', trimmed, + '--readFilesCommand', 'zcat', '--outFilterMultimapNmax', '10', '--outFilterMismatchNoverLmax', '0.05', + '--outFilterScoreMinOverLread', '0.75', '--outFilterMatchNminOverLread', '0.85', '--alignIntronMax', '1', + '--outFilterIntronMotifs', 'RemoveNoncanonical', '--outSAMtype', 'BAM', 'SortedByCoordinate', + '--quantMode', 'GeneCounts', '--limitBAMsortRAM', '4000000000', '--genomeLoad', 'NoSharedMemory'], cwd=star) + bam = analysis / 'BAM' / (sample + '.bam') + shutil.move(star / 'Aligned.sortedByCoord.out.bam', bam) + run(['samtools', 'quickcheck', bam]) + run(['samtools', 'index', '-@', args.threads, bam]) + metrics = dict(line.strip().split('|', 1) for line in (star / 'Log.final.out').read_text().splitlines() if '|' in line) + metrics = {k.strip(): v.strip() for k, v in metrics.items()} + logs.append([sample] + [metrics[k] for k in ['Number of input reads', 'Uniquely mapped reads number', 'Number of reads mapped to multiple loci', 'Number of reads mapped to too many loci']]) + for feature, mode in [('exon', 'intersection-strict'), ('intron', 'union')]: + annotation = ref / ('hg38dm6sacCer3.' + ('exons' if feature == 'exon' else 'introns') + '.gtf') + run(['htseq-count', '--format=bam', '--order=pos', '--stranded=reverse', '--minaqual=10', '--type=exon', + '--idattr=gene_id', '--mode=' + mode, '--quiet', bam, annotation], + output=analysis / 'htseq_exons_introns' / f'{sample}_{feature}.tsv') + (analysis / 'Logs/Logs.tsv').write_text(''.join(' '.join(row) + '\n' for row in logs)) + run(['perl', HOME / 'Scripts/htseqGeneCountsMerge.pl'], cwd=analysis / 'htseq_exons_introns') + env = os.environ.copy() + env.update(GTF_NO_EXON=str(ref / 'hg38dm6sacCer3.exons.no.gtf'), GTF_INTRON=str(ref / 'hg38dm6sacCer3.introns.gtf'), + GENOME_SIZES=str(ref / 'sizes.genome'), BAM_PATH=str(analysis / 'BAM')) + run(['bash', HOME / 'preDRUID.sh'], cwd=analysis, env=env) + shutil.copyfile(args.manifest, analysis / 'samples.tsv') + shutil.copyfile('/opt/debian-packages.tsv', analysis / 'software-versions.tsv') + shutil.copyfile(ref / 'provenance.json', analysis / 'reference-provenance.json') + (analysis / 'workflow-provenance.json').write_text(json.dumps({ + 'upstream_commit': (HOME / 'UPSTREAM_COMMIT').read_text().strip(), + 'script_sha256': {p.name: digest(p) for p in [HERE / 'druid.py', HERE / 'fit.R', + HOME / 'DRUID.R', HOME / 'coverage.py', HOME / 'preDRUID.sh', HOME / 'process_GTF.R']}, + }, indent=2) + '\n') + fit(args) + + +def fit(args): + analysis = args.work / 'Analysis' / args.dataset + run(['Rscript', HERE / 'fit.R', analysis, args.manifest]) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('command', choices=['download', 'reference', 'analyze', 'fit', 'all', 'smoke-test']) + parser.add_argument('--work', type=Path, default=Path.cwd(), help='Mounted working directory (default: current directory)') + parser.add_argument('--manifest', type=Path, default=Path('/opt/config/HEK293_R1.tsv')) + parser.add_argument('--dataset', default='HEK293_R1') + parser.add_argument('--threads', type=int, default=int(os.environ.get('THREADS', '8'))) + parser.add_argument('--fasta', type=Path, help='Already species-labeled FASTA, with --gtf') + parser.add_argument('--gtf', type=Path, help='Already species-labeled GTF, with --fasta') + parser.add_argument('--small-genome', action='store_true', help='Small synthetic reference only; never use for the full human genome') + args = parser.parse_args() + args.work = args.work.resolve() + args.manifest = args.manifest.resolve() + if args.fasta: args.fasta = args.fasta.resolve() + if args.gtf: args.gtf = args.gtf.resolve() + if args.threads < 1 or not re.fullmatch(r'[A-Za-z0-9_.-]+', args.dataset) or args.dataset in ('.', '..'): + parser.error('Threads must be positive and dataset must be a simple directory name') + args.work.mkdir(parents=True, exist_ok=True) + Path(os.environ.get('TMPDIR', str(args.work / 'tmp'))).mkdir(parents=True, exist_ok=True) + if args.command == 'all': + memory_check() + fetch_fastq(args) + prepare_reference(args) + analyze(args) + elif args.command == 'smoke-test': + run(['bash', '/opt/tests/smoke-test.sh', args.work]) + else: + {'download': fetch_fastq, 'reference': prepare_reference, 'analyze': analyze, 'fit': fit}[args.command](args) + + +if __name__ == '__main__': + try: + main() + except (ValueError, RuntimeError, OSError, subprocess.CalledProcessError) as error: + sys.exit(f'ERROR: {error}') diff --git a/images/druid/complete/scripts/fit.R b/images/druid/complete/scripts/fit.R new file mode 100644 index 0000000..bb2fa27 --- /dev/null +++ b/images/druid/complete/scripts/fit.R @@ -0,0 +1,49 @@ +#!/usr/bin/env Rscript +args <- commandArgs(trailingOnly = TRUE) +stopifnot(length(args) == 2) +source(file.path(Sys.getenv('DRUID_HOME', '/opt/DRUID'), 'DRUID.R')) +samples <- read.delim(args[2], check.names = FALSE, stringsAsFactors = FALSE) +samples <- samples[order(samples$time_hours), ] +stopifnot(!anyDuplicated(samples$sample), !anyDuplicated(samples$time_hours)) +setwd(args[1]) +times <- samples$time_hours +htseq <- read.delim('htseq_exons_introns/combined_exon.txt', check.names = FALSE) +stopifnot(all(samples$sample %in% names(htseq))) +get_counts <- function(organism) { + selected <- htseq$Organism == organism + counts <- htseq[selected, samples$sample, drop = FALSE] + rownames(counts) <- htseq$Gene[selected] + Filter(counts) +} +data.reads <- get_counts('h') +spike.reads <- get_counts('f') +if (nrow(data.reads) < 2 || nrow(spike.reads) < 1) stop('Insufficient human or fly genes passing the tutorial expression filter') +logs <- read.table('Logs/Logs.tsv', col.names = c('Sample','Input','Unique','Multi','TooMany')) +stopifnot(!anyDuplicated(logs$Sample)) +library.sizes <- logs$Unique[match(samples$sample, logs$Sample)] +if (anyNA(library.sizes) || any(library.sizes <= 0)) stop('Missing or zero uniquely mapped read counts') +# Explicit sample/type order avoids lexical time-point and filesystem-order errors. +csvs <- unlist(lapply(samples$sample, function(sample) { + file.path('intersect', paste0(sample, c('_exon_neg.csv','_exon_pos.csv','_intron_neg.csv','_intron_pos.csv'))) +})) +stopifnot(all(file.exists(csvs))) +suffixes <- paste0('.', times) +exons.all <- rbind(ExtractAndMerge(csvs, 4, 0, suffixes), ExtractAndMerge(csvs, 4, 1, suffixes)) +introns.all <- rbind(ExtractAndMerge(csvs, 4, 2, suffixes), ExtractAndMerge(csvs, 4, 3, suffixes)) +introns.human <- introns.all[introns.all$organism == 'h', ] +coverage <- as.matrix(introns.human[, grep('^mean\\.', names(introns.human)), drop=FALSE]) +eligible <- complete.cases(coverage) & apply(coverage, 1, function(x) all(x >= 0.5)) +if (sum(eligible) < 5) stop('Fewer than five introns pass the coverage threshold; more reads or suitable introns are required') +spike.norm <- colSums(spike.reads) +stopifnot(all(is.finite(spike.norm)), all(spike.norm > 0)) +half.lives.spikeins <- LeaveOneOut(times, data.reads, spike.norm, HalfLives4SU) +half.lives.DRUID <- DRUID(x=times, introns=introns.human, data=data.reads, + lib.sizes=library.sizes, path='DRUID.replicate.1', kmeans.seed=42) +write.csv(half.lives.spikeins, 'half-lives.spikeins.csv') +write.csv(half.lives.DRUID, 'half-lives.DRUID.csv') +saveRDS(list(samples=samples, library.sizes=library.sizes, data.reads=data.reads, + spike.reads=spike.reads, exons=exons.all, introns=introns.all, + half.lives.DRUID=half.lives.DRUID, half.lives.spikeins=half.lives.spikeins), + 'analysis.rds') +writeLines(capture.output(sessionInfo()), 'R-sessionInfo.txt') +cat('Saved half-lives, heatmap, intron normalization counts, and analysis.rds\n') diff --git a/images/druid/complete/scripts/patch_upstream.py b/images/druid/complete/scripts/patch_upstream.py new file mode 100644 index 0000000..d006f16 --- /dev/null +++ b/images/druid/complete/scripts/patch_upstream.py @@ -0,0 +1,48 @@ +"""Small, checked patches to pinned upstream; the original stays in *.upstream.*.""" +from pathlib import Path +import shutil +root = Path('/opt/DRUID') +p = root / 'DRUID.R' +shutil.copy(p, root / 'DRUID.upstream.R') +s = p.read_text() +s = s[:s.index('# Example')] +old = ''' introns.mean.filt.norm <- introns.mean.filt / lib.sizes + maxes <- apply(introns.mean.filt.norm, 1, max) + maxes.matrix <- matrix(rep(maxes, length(maxes)), ncol = length(maxes)) + introns.mean.filt.norm <- introns.mean.filt.norm / maxes.matrix''' +assert s.count(old) == 1 +s = s.replace(old, ''' introns.mean.filt.norm <- sweep(as.matrix(introns.mean.filt), 2, lib.sizes, "/") + maxes <- apply(introns.mean.filt.norm, 1, max) + introns.mean.filt.norm <- sweep(introns.mean.filt.norm, 1, maxes, "/")''') +p.write_text(s) +# coverage.py was Python 2 with mixed tabs; preserve numerical coverage semantics. +p = root / 'coverage.py' +shutil.copy(p, root / 'coverage.upstream.py') +shutil.copy('/opt/workflow/coverage.py', p) +# Replace the hard-coded-path shell wrapper with a fail-fast equivalent. +p = root / 'preDRUID.sh' +shutil.copy(p, root / 'preDRUID.upstream.sh') +shutil.copy('/opt/workflow/preDRUID.sh', p) +# GenomicFeatures renamed the exons() selection argument from vals to filter. +p = root / 'process_GTF.R' +shutil.copy(p, root / 'process_GTF.upstream.R') +s = p.read_text() +old = 'vals = list(tx_id=transcripts.length$tx_id)' +assert s.count(old) == 1 +p.write_text(s.replace(old, 'filter = list(tx_id=transcripts.length$tx_id)')) +# Modern GenomicRanges no longer accepts ignoreSelf in this method. Remove +# diagonal self hits explicitly, preserving strand-aware exon overlap filtering. +p = root / 'process_GTF.R' +s = p.read_text() +old = '''ov <- findOverlaps( + exons, + type = "any", + ignoreSelf = TRUE, + select = "arbitrary" +) +exons <- exons[is.na(ov)]''' +assert s.count(old) == 1 +s = s.replace(old, '''ov <- findOverlaps(exons, exons, type = "any") +overlapping <- queryHits(ov)[queryHits(ov) != subjectHits(ov)] +exons <- exons[!seq_along(exons) %in% overlapping]''') +p.write_text(s) diff --git a/images/druid/complete/scripts/preDRUID.sh b/images/druid/complete/scripts/preDRUID.sh new file mode 100644 index 0000000..30685bc --- /dev/null +++ b/images/druid/complete/scripts/preDRUID.sh @@ -0,0 +1,31 @@ +#!/usr/bin/env bash +set -euo pipefail +: "${GTF_NO_EXON:?}" "${GTF_INTRON:?}" "${GENOME_SIZES:?}" "${BAM_PATH:?}" +export LC_ALL=C +mkdir -p "${TMPDIR:-/tmp}" +scratch=$(mktemp -d "${TMPDIR:-/tmp}/predruid.XXXXXX") +trap 'rm -rf "$scratch"' EXIT +for type in exon intron; do + input=$GTF_NO_EXON + [[ $type == intron ]] && input=$GTF_INTRON + awk '$7=="-"' "$input" | sort -k1,1 -k4,4n > "$scratch/$type.neg.gtf" + awk '$7=="+"' "$input" | sort -k1,1 -k4,4n > "$scratch/$type.pos.gtf" +done +mkdir -p intersect +shopt -s nullglob +bams=("$BAM_PATH"/*.bam) +(( ${#bams[@]} > 0 )) || { echo 'No BAM files found' >&2; exit 1; } +for bam in "${bams[@]}"; do + base=$(basename "$bam" .bam) + samtools view -h "$bam" | awk '/^@/ || /\tNH:i:1(\t|$)/' | samtools view -b -o "$scratch/unique.bam" - + bedtools genomecov -ibam "$scratch/unique.bam" -bg -split -strand - | sort -k1,1 -k2,2n > "$scratch/neg.bg" + bedtools genomecov -ibam "$scratch/unique.bam" -bg -split -strand + | sort -k1,1 -k2,2n > "$scratch/pos.bg" + for type in exon intron; do + for strand in neg pos; do + other=pos + [[ $strand == pos ]] && other=neg + bedtools intersect -a "$scratch/$type.$other.gtf" -b "$scratch/$strand.bg" -wo -sorted -g "$GENOME_SIZES" > "$scratch/overlap.tsv" + python3 "${DRUID_HOME:-/opt/DRUID}/coverage.py" "$scratch/overlap.tsv" > "intersect/${base}_${type}_${strand}.csv" + done + done +done diff --git a/images/druid/complete/tests/check_annotations.py b/images/druid/complete/tests/check_annotations.py new file mode 100644 index 0000000..b170497 --- /dev/null +++ b/images/druid/complete/tests/check_annotations.py @@ -0,0 +1,29 @@ +#!/usr/bin/env python3 +"""Biological invariants of representative-transcript and overlap filtering.""" +from pathlib import Path +import subprocess +import tempfile + +with tempfile.TemporaryDirectory() as temp: + root = Path(temp) + gtf = root / 'input.gtf' + rows = [] + for gene, tx, strand, intervals in [ + ('A', 'AS', '+', [(101,150),(301,350)]), + ('A', 'AL', '+', [(101,200),(301,400)]), + ('B', 'B', '+', [(181,220),(501,550)]), + ('C', 'C', '-', [(101,200),(301,400)]), + ]: + for i,(a,b) in enumerate(intervals): + n = i+1 if strand=='+' else 2-i + attrs = f'gene_id "{gene}h"; transcript_id "{tx}h"; exon_number {n}; exon_id "{tx}.{n}h"; gene_name "{gene}h";' + rows.append(f'chr1h\trefGene\texon\t{a}\t{b}\t.\t{strand}\t.\t{attrs}\n') + gtf.write_text(''.join(rows)) + subprocess.run(['Rscript','/opt/DRUID/process_GTF.R',str(gtf)], check=True) + def features(suffix): + rows = [x.split('\t') for x in (root / ('input.'+suffix+'.gtf')).read_text().splitlines() if not x.startswith('#')] + return {(int(r[3]),int(r[4]),r[6]) for r in rows} + assert features('exons') == {(101,200,'+'),(301,400,'+'),(181,220,'+'),(501,550,'+'),(101,200,'-'),(301,400,'-')} + assert features('exons.no') == {(301,400,'+'),(501,550,'+'),(101,200,'-'),(301,400,'-')} + assert features('introns') == {(201,300,'-')} +print('PASS: longest transcript, exon overlaps, intron exclusion, and strand handling') diff --git a/images/druid/complete/tests/check_coverage.py b/images/druid/complete/tests/check_coverage.py new file mode 100644 index 0000000..5ac093f --- /dev/null +++ b/images/druid/complete/tests/check_coverage.py @@ -0,0 +1,19 @@ +#!/usr/bin/env python3 +"""Test overlap coordinates, zero coverage, and feature metadata against hand calculations.""" +import csv +import io +import subprocess +import tempfile +import numpy as np +attrs = 'gene_id "Ah"; transcript_id "TXh"; exon_number "1"; exon_id "TX.1h"; gene_name "Ah";' +# Four-base feature; coverage [2, 2, 3, 0]. Input intervals use BED coordinates. +rows = [f'chr1h\tx\texon\t11\t14\t.\t+\t.\t{attrs}\tchr1h\t10\t12\t2\t2', + f'chr1h\tx\texon\t11\t14\t.\t+\t.\t{attrs}\tchr1h\t12\t13\t3\t1'] +with tempfile.NamedTemporaryFile(mode='w') as f: + f.write('\n'.join(rows)+'\n'); f.flush() + output = subprocess.check_output(['python3','/opt/DRUID/coverage.py',f.name], text=True) +row = list(csv.DictReader(io.StringIO(output)))[0] +assert row['gene']=='A' and row['transcript']=='TX' and row['size']=='4' +assert float(row['mean'])==1.75 and float(row['coverage'])==0.75 +assert np.isclose(float(row['stdev']), np.std([2,2,3,0])) +print('PASS: coverage coordinates and statistics') diff --git a/images/druid/complete/tests/check_predruid.py b/images/druid/complete/tests/check_predruid.py new file mode 100644 index 0000000..38679e2 --- /dev/null +++ b/images/druid/complete/tests/check_predruid.py @@ -0,0 +1,27 @@ +#!/usr/bin/env python3 +"""Ensure preDRUID excludes NH=10 reads and selects the opposite read strand.""" +import csv +import os +from pathlib import Path +import subprocess +import tempfile +with tempfile.TemporaryDirectory() as temp: + root = Path(temp) + bamdir = root/'BAM'; bamdir.mkdir() + attrs = 'gene_id "Ah"; transcript_id "TXh"; exon_number 1; exon_id "TX.1h"; gene_name "Ah";' + (root/'exons.gtf').write_text(f'chr1h\tx\texon\t11\t14\t.\t+\t.\t{attrs}\n') + (root/'introns.gtf').write_text(f'chr1h\tx\texon\t21\t24\t.\t+\t.\t{attrs}\n') + (root/'sizes').write_text('chr1h\t100\n') + sam = '@HD\tVN:1.6\tSO:coordinate\n@SQ\tSN:chr1h\tLN:100\n' + for i,(flag,pos,nh) in enumerate([(16,11,1),(16,11,10),(0,11,1),(16,21,1),(16,21,10)]): + sam += f'r{i}\t{flag}\tchr1h\t{pos}\t255\t4M\t*\t0\t0\tAAAA\tIIII\tNH:i:{nh}\n' + (root/'test.sam').write_text(sam) + subprocess.run(['samtools','view','-b','-o',str(bamdir/'test.bam'),str(root/'test.sam')],check=True) + env = dict(os.environ,GTF_NO_EXON=str(root/'exons.gtf'),GTF_INTRON=str(root/'introns.gtf'), + GENOME_SIZES=str(root/'sizes'),BAM_PATH=str(bamdir),TMPDIR=temp) + subprocess.run(['bash','/opt/DRUID/preDRUID.sh'],cwd=root,env=env,check=True) + for feature in ['exon','intron']: + rows=list(csv.DictReader((root/f'intersect/test_{feature}_neg.csv').open())) + assert len(rows)==1 and float(rows[0]['mean'])==1 and float(rows[0]['coverage'])==1 + assert not list(csv.DictReader((root/f'intersect/test_{feature}_pos.csv').open())) +print('PASS: exact NH=1 filtering and reverse-stranded coverage') diff --git a/images/druid/complete/tests/check_results.R b/images/druid/complete/tests/check_results.R new file mode 100644 index 0000000..35041b1 --- /dev/null +++ b/images/druid/complete/tests/check_results.R @@ -0,0 +1,21 @@ +args <- commandArgs(trailingOnly=TRUE) +a <- readRDS(file.path(args[1], 'Analysis/HEK293_R1/analysis.rds')) +stopifnot(ncol(a$half.lives.DRUID)==7, nrow(a$half.lives.DRUID)==48, + sum(is.finite(a$half.lives.DRUID[,1])) > 30, + nrow(a$half.lives.spikeins)==48, + all(a$library.sizes > 0), nrow(a$introns)==60) +source('/opt/DRUID/DRUID.R') +x <- c(1,2,4,8,12,24) +y <- 100 * (1-exp(-log(2)*x/3)) +stopifnot(abs(HalfLives4SU(x,y,doubling.time=0)-3) < 1e-4) +# Changing sequencing depth per sample must preserve clustering and half-lives. +introns <- a$introns[a$introns$organism=='h',] +original <- DRUID(x, introns, a$data.reads, a$library.sizes, + tempfile(), plot=FALSE, kmeans.seed=42) +factors <- c(2,5,3,4,7,6) +mi <- grep('^mean\\.',names(introns)) +introns[,mi] <- sweep(introns[,mi],2,factors,'*') +scaled <- DRUID(x, introns, sweep(a$data.reads,2,factors,'*'), a$library.sizes*factors, + tempfile(), plot=FALSE, kmeans.seed=42) +stopifnot(isTRUE(all.equal(original,scaled,tolerance=1e-5))) +cat('PASS: end-to-end outputs, known half-life, and sequencing-depth invariance\n') diff --git a/images/druid/complete/tests/make_fixture.py b/images/druid/complete/tests/make_fixture.py new file mode 100644 index 0000000..efdfb02 --- /dev/null +++ b/images/druid/complete/tests/make_fixture.py @@ -0,0 +1,55 @@ +#!/usr/bin/env python3 +"""Deterministic stranded, six-time-point reads for end-to-end execution testing.""" +import csv +import gzip +import math +from pathlib import Path +import random +import sys + +root = Path(sys.argv[1]) +root.mkdir(parents=True, exist_ok=True) +rng = random.Random(99517) +sequences = {name: ''.join(rng.choices('ACGT', k=length)) for name, length in [('chr1h', 60000), ('chr2Lf', 8000), ('chrIy', 8000)]} +with (root / 'fixture.fa').open('w') as handle: + for name, seq in sequences.items(): + handle.write(f'>{name}\n{seq}\n') +genes = [] +with (root / 'fixture.gtf').open('w') as handle: + for chrom, number in [('chr1h', 48), ('chr2Lf', 6), ('chrIy', 6)]: + suffix = chrom[-1] + for i in range(number): + start = i * 1000 + 100 + strand = '+' if i % 2 == 0 else '-' + gene, transcript = f'GENE{i:03}{suffix}', f'TX{i:03}{suffix}' + genes.append((chrom, i, start, strand)) + for j, (a, b) in enumerate([(start, start+199), (start+400, start+599)]): + exon = j+1 if strand == '+' else 2-j + attrs = f'gene_id "{gene}"; transcript_id "{transcript}"; exon_number {exon}; exon_id "TX{i:03}.{exon}{suffix}"; gene_name "{gene}";' + handle.write(f'{chrom}\trefGene\texon\t{a}\t{b}\t.\t{strand}\t.\t{attrs}\n') +raw = root / 'RawSequencingData/HEK293_R1' +raw.mkdir(parents=True, exist_ok=True) +times = [1,2,4,8,12,24] +with (root / 'samples.tsv').open('w') as handle: + writer = csv.writer(handle, delimiter='\t') + writer.writerow(['sample','time_hours']) + for t in times: + sample = f'HEK293_R1_{t:02}' + writer.writerow([sample,t]) + counter = 0 + with gzip.open(raw / (sample+'.fastq.gz'), 'wt') as fastq: + for chrom, i, start, strand in genes: + human = chrom.endswith('h') + exonic = round(180 * (1 - math.exp(-math.log(2) * t / (1+i%12)))) if human else 80 + # Four distinct intron dynamics, with a well-behaved early plateau. + shape = [1, 0.05 + t/24, 0.05 + (t/24)**2, 0.05 + (t/24)**3][i%4] + intronic = round((65+i) * shape) if human else 20 + for kind, nreads in [('exon',exonic), ('intron',intronic)]: + for k in range(nreads): + pos = start-1 + (0 if kind=='exon' else 200) + ((k*13) % 151) + read = sequences[chrom][pos:pos+50] + if strand == '+': + read = read.translate(str.maketrans('ACGT','TGCA'))[::-1] + counter += 1 + fastq.write(f'@{sample}_{counter}\n{read}\n+\n'+ 'I'*50+'\n') +print('Synthetic fixture created:', root) diff --git a/images/druid/complete/tests/smoke-test.sh b/images/druid/complete/tests/smoke-test.sh new file mode 100644 index 0000000..ef4b776 --- /dev/null +++ b/images/druid/complete/tests/smoke-test.sh @@ -0,0 +1,41 @@ +#!/usr/bin/env bash +set -euo pipefail +root=$(realpath "${1:-/work}") +mkdir -p "$root" +fixture=$(mktemp -d "$root/smoke.XXXXXX") +export TMPDIR="$fixture/tmp" +mkdir -p "$TMPDIR" +python3 /opt/tests/check_coverage.py +python3 /opt/tests/check_predruid.py +python3 /opt/tests/check_annotations.py +python3 /opt/tests/make_fixture.py "$fixture" +druid reference --work "$fixture" --fasta "$fixture/fixture.fa" --gtf "$fixture/fixture.gtf" --small-genome --threads 2 +druid analyze --work "$fixture" --manifest "$fixture/samples.tsv" --threads 2 +Rscript /opt/tests/check_results.R "$fixture" +test -s "$fixture/Analysis/HEK293_R1/DRUID.replicate.1.heatmap.DRUID.pdf" +# Also exercise UCSC converters using a small real genePred table and 2bit genome. +curl -fL --retry 3 https://hgdownload.soe.ucsc.edu/goldenPath/sacCer3/database/sgdGene.txt.gz -o "$fixture/sgdGene.txt.gz" +python3 - "$fixture" <<'PYCODE' +import sys +from pathlib import Path +sys.path.insert(0, '/opt/workflow') +from druid import unpack_genepred +root = Path(sys.argv[1]) +unpack_genepred(root / 'sgdGene.txt.gz', root / 'sgdGene.gp', 'sgdGene') +PYCODE +genePredToGtf -utr file "$fixture/sgdGene.gp" "$fixture/sgdGene.gtf" +test -s "$fixture/sgdGene.gtf" +curl -fL --retry 3 https://hgdownload.soe.ucsc.edu/goldenPath/sacCer3/bigZips/sacCer3.2bit -o "$fixture/sacCer3.2bit" +twoBitToFa "$fixture/sacCer3.2bit" "$fixture/sacCer3.fa" +test -s "$fixture/sacCer3.fa" +python3 - "$fixture" <<'PYCODE' +import sys +from pathlib import Path +sys.path.insert(0, '/opt/workflow') +from druid import label_reference +root = Path(sys.argv[1]) +label_reference(root / 'sacCer3.fa', root / 'sgdGene.gtf', 'y', False, root / 'yeast.fa', root / 'yeast.gtf') +PYCODE +Rscript /opt/DRUID/process_GTF.R "$fixture/yeast.gtf" +test -s "$fixture/yeast.introns.gtf" +echo "PASS: complete workflow and UCSC converters. Results: $fixture"