Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
34 lines
1.9 KiB
Bash
34 lines
1.9 KiB
Bash
#!/bin/bash
|
|
# v2 methylation QC: numerics vs the sesame oracle (R, MIT) on the public minfiDataEPIC set (3 x NA12878, EPIC v1).
|
|
# Args: <unused data dir> <results dir>. Run with OMNI=0.
|
|
set -uo pipefail
|
|
OUT=$2
|
|
R=$OUT/meth; W=/work/meth; mkdir -p $R $W
|
|
ARCH=$(uname -m)
|
|
echo "=== data $(date -u +%T) ==="
|
|
cd $W
|
|
curl -sSfL -O https://bioconductor.org/packages/release/data/experiment/src/contrib/minfiDataEPIC_1.38.0.tar.gz \
|
|
&& tar -xzf minfiDataEPIC_1.38.0.tar.gz minfiDataEPIC/inst/extdata
|
|
curl -sSfL -O https://webdata.illumina.com/downloads/productfiles/methylationEPIC/infinium-methylationepic-v-1-0-b5-manifest-file-csv.zip \
|
|
&& unzip -oq infinium-methylationepic-v-1-0-b5-manifest-file-csv.zip
|
|
cd - > /dev/null
|
|
MAN=$W/infinium-methylationepic-v-1-0-b5-manifest-file.csv
|
|
E=$W/minfiDataEPIC/inst/extdata
|
|
PREFIXES="$E/200144450018/200144450018_R04C01 $E/200144450019/200144450019_R07C01 $E/200144450021/200144450021_R05C01"
|
|
.venv/bin/simurg methyl convert-manifest --manifest $MAN --out $W/mft --quiet
|
|
|
|
echo "=== sesame (conda) $(date -u +%T) ==="
|
|
curl -sSfL -o $W/mf.sh "https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-$ARCH.sh" \
|
|
&& bash $W/mf.sh -b -p $W/conda > /dev/null
|
|
$W/conda/bin/mamba create -y -q -p $W/renv -c conda-forge -c bioconda r-base bioconductor-sesame > $R/conda.log 2>&1 \
|
|
|| { echo "!! sesame install failed"; tail -30 $R/conda.log; }
|
|
$W/renv/bin/Rscript scripts/v2/sesame_oracle.R $W/mft/probes.tsv $W/mft/controls.tsv $W/sesame $PREFIXES \
|
|
> $R/sesame.log 2>&1 || echo "!! sesame run failed"
|
|
tail -40 $R/sesame.log
|
|
|
|
echo "=== compare $(date -u +%T) ==="
|
|
.venv/bin/python scripts/v2/methyl_oracle_compare.py --manifest $MAN --sesame $W/sesame --out $R/compare.json $PREFIXES
|
|
.venv/bin/simurg methyl qc --manifest $MAN --idats $E --sheet $E/SampleSheet.csv --out $W/qc --threads 0 --quiet
|
|
cp $W/qc/methyl_qc_summary.tsv $W/qc/controls_qc.tsv $R/
|
|
echo "=== task done $(date -u +%T) ==="
|