Add self-destructing AWS trial and large-scale P1 validation script
scripts/cloud/trial.sh launches one spot instance that validates the parsers on the 1000 Genomes Omni2.5 set (2141 samples, vendor GTCs) streamed from the public S3 bucket, uploads the summary and terminates itself. Cost guards: spot price cap, hard time limit (shutdown = terminate), EXIT-trap shutdown, CloudWatch idle alarm that terminates from outside, disk deleted with the instance, monthly budget alert, and an `audit` command listing anything still billable. scripts/validate_p1_large.py compares every GTC with BeadArrayFiles, checks byte-identical round-trips, IDAT pairs and BPM/EGT/CSV. Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
This commit is contained in:
1 parent
4357d3deef
commit
a9aa6dc080
4 files changed
+412
No files matched your search
@@ -16,3 +16,5 @@ __pycache__/
|
|||||||
.ruff_cache/
|
.ruff_cache/
|
||||||
# rust
|
# rust
|
||||||
target/
|
target/
|
||||||
|
.trial-last
|
||||||
|
results/
|
||||||
Executable
+123
@@ -0,0 +1,123 @@
|
|||||||
|
#!/usr/bin/env bash
|
||||||
|
# Self-destructing, budget-guarded AWS validation run for SimurgArray.
|
||||||
|
#
|
||||||
|
# bash scripts/cloud/trial.sh setup <alert-email> # one-time: budget alert, results bucket, instance role
|
||||||
|
# bash scripts/cloud/trial.sh launch # start one spot run (terminates itself when done)
|
||||||
|
# bash scripts/cloud/trial.sh status # instance state + last lines of the remote log
|
||||||
|
# bash scripts/cloud/trial.sh results # download log/summary of the latest run
|
||||||
|
# bash scripts/cloud/trial.sh kill # emergency: terminate every trial instance now
|
||||||
|
# bash scripts/cloud/trial.sh audit # list anything still costing money
|
||||||
|
#
|
||||||
|
# Cost guards: spot price cap, hard time cap (shutdown = terminate), EXIT trap shutdown, CloudWatch idle alarm
|
||||||
|
# that terminates from outside, disk deleted with the instance, monthly budget e-mail alert.
|
||||||
|
set -euo pipefail
|
||||||
|
cd "$(dirname "$0")/../.."
|
||||||
|
REGION=${REGION:-us-east-1}
|
||||||
|
TYPE=${TYPE:-m7g.4xlarge}
|
||||||
|
MAX_PRICE=${MAX_PRICE:-0.30} # USD/hour spot cap (on-demand is 0.6528)
|
||||||
|
MAX_MINUTES=${MAX_MINUTES:-360} # hard cap: 6 hours
|
||||||
|
DISK_GB=${DISK_GB:-200}
|
||||||
|
BUDGET_USD=${BUDGET_USD:-10}
|
||||||
|
TAG=simurg-trial
|
||||||
|
ACCOUNT=$(aws sts get-caller-identity --query Account --output text)
|
||||||
|
BUCKET=simurg-trial-$ACCOUNT
|
||||||
|
ROLE=simurg-trial-role
|
||||||
|
AWS="aws --region $REGION"
|
||||||
|
|
||||||
|
cmd=${1:-help}
|
||||||
|
case "$cmd" in
|
||||||
|
setup)
|
||||||
|
EMAIL=${2:?usage: trial.sh setup <alert-email>}
|
||||||
|
echo "-> budget alert: \$$BUDGET_USD/month, e-mail $EMAIL at 50%, 100% actual and 100% forecast"
|
||||||
|
cat > /tmp/simurg-budget.json <<EOF
|
||||||
|
{"BudgetName":"$TAG-monthly","BudgetLimit":{"Amount":"$BUDGET_USD","Unit":"USD"},"TimeUnit":"MONTHLY","BudgetType":"COST"}
|
||||||
|
EOF
|
||||||
|
notif() { echo "{\"Notification\":{\"NotificationType\":\"$1\",\"ComparisonOperator\":\"GREATER_THAN\",\"Threshold\":$2,\"ThresholdType\":\"PERCENTAGE\"},\"Subscribers\":[{\"SubscriptionType\":\"EMAIL\",\"Address\":\"$EMAIL\"}]}"; }
|
||||||
|
echo "[$(notif ACTUAL 50),$(notif ACTUAL 100),$(notif FORECASTED 100)]" > /tmp/simurg-notif.json
|
||||||
|
aws budgets create-budget --account-id "$ACCOUNT" --budget file:///tmp/simurg-budget.json \
|
||||||
|
--notifications-with-subscribers file:///tmp/simurg-notif.json 2>/dev/null || echo " (budget already exists)"
|
||||||
|
echo "-> results bucket s3://$BUCKET (objects expire after 30 days)"
|
||||||
|
$AWS s3api create-bucket --bucket "$BUCKET" >/dev/null 2>&1 || true
|
||||||
|
$AWS s3api put-public-access-block --bucket "$BUCKET" --public-access-block-configuration \
|
||||||
|
BlockPublicAcls=true,IgnorePublicAcls=true,BlockPublicPolicy=true,RestrictPublicBuckets=true
|
||||||
|
$AWS s3api put-bucket-lifecycle-configuration --bucket "$BUCKET" --lifecycle-configuration \
|
||||||
|
'{"Rules":[{"ID":"expire","Status":"Enabled","Filter":{},"Expiration":{"Days":30}}]}'
|
||||||
|
echo "-> instance role $ROLE (read/write only this bucket)"
|
||||||
|
aws iam create-role --role-name $ROLE --assume-role-policy-document \
|
||||||
|
'{"Version":"2012-10-17","Statement":[{"Effect":"Allow","Principal":{"Service":"ec2.amazonaws.com"},"Action":"sts:AssumeRole"}]}' \
|
||||||
|
>/dev/null 2>&1 || true
|
||||||
|
aws iam put-role-policy --role-name $ROLE --policy-name bucket-only --policy-document \
|
||||||
|
"{\"Version\":\"2012-10-17\",\"Statement\":[{\"Effect\":\"Allow\",\"Action\":[\"s3:GetObject\",\"s3:PutObject\"],\"Resource\":\"arn:aws:s3:::$BUCKET/*\"}]}"
|
||||||
|
aws iam create-instance-profile --instance-profile-name $ROLE >/dev/null 2>&1 || true
|
||||||
|
aws iam add-role-to-instance-profile --instance-profile-name $ROLE --role-name $ROLE 2>/dev/null || true
|
||||||
|
echo "setup done (budget alerts need you to confirm nothing; AWS e-mails $EMAIL when thresholds are crossed)"
|
||||||
|
;;
|
||||||
|
launch)
|
||||||
|
running=$($AWS ec2 describe-instances --filters Name=tag:Project,Values=$TAG \
|
||||||
|
Name=instance-state-name,Values=pending,running --query 'length(Reservations[].Instances[])')
|
||||||
|
[[ "$running" == "0" ]] || { echo "a trial instance is already running; refusing to start another"; exit 1; }
|
||||||
|
RUN=run-$(date -u +%Y%m%dT%H%M%SZ)
|
||||||
|
echo "-> packing repository at $(git rev-parse --short HEAD) -> s3://$BUCKET/$RUN/"
|
||||||
|
git archive --format=tar.gz -o /tmp/simurg-repo.tar.gz HEAD
|
||||||
|
$AWS s3 cp /tmp/simurg-repo.tar.gz "s3://$BUCKET/$RUN/repo.tar.gz" --only-show-errors
|
||||||
|
sed -e "s/__BUCKET__/$BUCKET/" -e "s/__RUN__/$RUN/" -e "s/__MAX_MINUTES__/$MAX_MINUTES/" \
|
||||||
|
scripts/cloud/userdata.sh > /tmp/simurg-userdata.sh
|
||||||
|
AMI=$($AWS ssm get-parameter --name /aws/service/canonical/ubuntu/server/24.04/stable/current/arm64/hvm/ebs-gp3/ami-id \
|
||||||
|
--query Parameter.Value --output text)
|
||||||
|
echo "-> launching $TYPE spot (cap \$$MAX_PRICE/h, hard limit ${MAX_MINUTES} min, ${DISK_GB} GB disk deleted on exit)"
|
||||||
|
ID=$($AWS ec2 run-instances --image-id "$AMI" --instance-type "$TYPE" --count 1 \
|
||||||
|
--iam-instance-profile Name=$ROLE \
|
||||||
|
--instance-initiated-shutdown-behavior terminate \
|
||||||
|
--instance-market-options "MarketType=spot,SpotOptions={MaxPrice=$MAX_PRICE,SpotInstanceType=one-time,InstanceInterruptionBehavior=terminate}" \
|
||||||
|
--block-device-mappings "DeviceName=/dev/sda1,Ebs={VolumeSize=$DISK_GB,VolumeType=gp3,DeleteOnTermination=true}" \
|
||||||
|
--metadata-options HttpTokens=required \
|
||||||
|
--user-data file:///tmp/simurg-userdata.sh \
|
||||||
|
--tag-specifications "ResourceType=instance,Tags=[{Key=Project,Value=$TAG},{Key=Run,Value=$RUN}]" \
|
||||||
|
"ResourceType=volume,Tags=[{Key=Project,Value=$TAG},{Key=Run,Value=$RUN}]" \
|
||||||
|
--query 'Instances[0].InstanceId' --output text)
|
||||||
|
echo "$ID $RUN" > .trial-last
|
||||||
|
echo "-> idle alarm: terminate if average CPU < 2% for 30 minutes"
|
||||||
|
$AWS cloudwatch put-metric-alarm --alarm-name "$TAG-idle-$ID" --namespace AWS/EC2 --metric-name CPUUtilization \
|
||||||
|
--dimensions Name=InstanceId,Value="$ID" --statistic Average --period 300 --evaluation-periods 6 \
|
||||||
|
--threshold 2 --comparison-operator LessThanThreshold --treat-missing-data notBreaching \
|
||||||
|
--alarm-actions "arn:aws:automate:$REGION:ec2:terminate"
|
||||||
|
echo "launched $ID ($RUN). Check with: bash scripts/cloud/trial.sh status"
|
||||||
|
;;
|
||||||
|
status)
|
||||||
|
read -r ID RUN < .trial-last
|
||||||
|
$AWS ec2 describe-instances --instance-ids "$ID" \
|
||||||
|
--query 'Reservations[0].Instances[0].[InstanceId,State.Name,InstanceType,LaunchTime]' --output text
|
||||||
|
$AWS s3 cp "s3://$BUCKET/$RUN/log.txt" - 2>/dev/null | tail -n "${LINES_:-15}" || echo "(no log uploaded yet)"
|
||||||
|
;;
|
||||||
|
results)
|
||||||
|
read -r ID RUN < .trial-last
|
||||||
|
mkdir -p "results/$RUN"
|
||||||
|
$AWS s3 cp "s3://$BUCKET/$RUN/" "results/$RUN/" --recursive --exclude repo.tar.gz --only-show-errors
|
||||||
|
ls -la "results/$RUN"
|
||||||
|
;;
|
||||||
|
kill)
|
||||||
|
ids=$($AWS ec2 describe-instances --filters Name=tag:Project,Values=$TAG \
|
||||||
|
Name=instance-state-name,Values=pending,running,stopping,stopped --query 'Reservations[].Instances[].InstanceId' --output text)
|
||||||
|
[[ -n "$ids" ]] && $AWS ec2 terminate-instances --instance-ids $ids --output text || echo "no trial instances"
|
||||||
|
;&
|
||||||
|
audit)
|
||||||
|
echo "-- instances (any project) --"
|
||||||
|
$AWS ec2 describe-instances --filters Name=instance-state-name,Values=pending,running,stopping,stopped \
|
||||||
|
--query 'Reservations[].Instances[].[InstanceId,State.Name,InstanceType]' --output text
|
||||||
|
echo "-- volumes --"
|
||||||
|
$AWS ec2 describe-volumes --query 'Volumes[].[VolumeId,State,Size]' --output text
|
||||||
|
echo "-- snapshots owned --"
|
||||||
|
$AWS ec2 describe-snapshots --owner-ids self --query 'Snapshots[].[SnapshotId,VolumeSize]' --output text
|
||||||
|
echo "-- elastic IPs --"
|
||||||
|
$AWS ec2 describe-addresses --query 'Addresses[].PublicIp' --output text
|
||||||
|
echo "-- trial alarms (deleting alarms of terminated instances) --"
|
||||||
|
for a in $($AWS cloudwatch describe-alarms --alarm-name-prefix "$TAG-idle-" --query 'MetricAlarms[].AlarmName' --output text); do
|
||||||
|
iid=${a#"$TAG-idle-"}
|
||||||
|
st=$($AWS ec2 describe-instances --instance-ids "$iid" --query 'Reservations[0].Instances[0].State.Name' --output text 2>/dev/null || echo gone)
|
||||||
|
if [[ "$st" == "terminated" || "$st" == "gone" ]]; then $AWS cloudwatch delete-alarms --alarm-names "$a"; echo "deleted $a"; else echo "$a ($st)"; fi
|
||||||
|
done
|
||||||
|
;;
|
||||||
|
*)
|
||||||
|
sed -n '2,15p' "$0"
|
||||||
|
;;
|
||||||
|
esac
|
||||||
Executable
+62
@@ -0,0 +1,62 @@
|
|||||||
|
#!/bin/bash
|
||||||
|
# EC2 user-data for a self-destructing SimurgArray validation run (Ubuntu 24.04, arm64).
|
||||||
|
# Placeholders (__X__) are filled in by scripts/cloud/trial.sh.
|
||||||
|
#
|
||||||
|
# Safety layers inside the instance:
|
||||||
|
# 1. `shutdown -h +MAX_MINUTES` is scheduled first -> hard time cap
|
||||||
|
# 2. an EXIT trap uploads the log and shuts down -> ends as soon as the job ends, success or failure
|
||||||
|
# Shutdown = TERMINATE (instance-initiated-shutdown-behavior) and the only disk has DeleteOnTermination=true,
|
||||||
|
# so nothing keeps costing money afterwards. A CloudWatch alarm (layer 3) terminates an idle instance from
|
||||||
|
# outside, in case the OS hangs.
|
||||||
|
set -uo pipefail
|
||||||
|
BUCKET="__BUCKET__"
|
||||||
|
RUN="__RUN__"
|
||||||
|
MAX_MINUTES="__MAX_MINUTES__"
|
||||||
|
shutdown -h +"$MAX_MINUTES" "SimurgArray trial hard time limit"
|
||||||
|
|
||||||
|
LOG=/var/log/simurg-trial.log
|
||||||
|
exec > >(tee -a "$LOG") 2>&1
|
||||||
|
S3=""
|
||||||
|
upload_log() { [[ -n "$S3" ]] && $S3 cp "$LOG" "s3://$BUCKET/$RUN/log.txt" --only-show-errors || true; }
|
||||||
|
finish() {
|
||||||
|
rc=$?
|
||||||
|
echo "=== finished rc=$rc at $(date -u +%FT%TZ) ==="
|
||||||
|
echo "$rc" > /tmp/rc && [[ -n "$S3" ]] && $S3 cp /tmp/rc "s3://$BUCKET/$RUN/exit_code.txt" --only-show-errors
|
||||||
|
upload_log
|
||||||
|
shutdown -h now
|
||||||
|
}
|
||||||
|
trap finish EXIT
|
||||||
|
set -e
|
||||||
|
echo "=== start $(date -u +%FT%TZ) on $(uname -m), $(nproc) cpus ==="
|
||||||
|
|
||||||
|
export DEBIAN_FRONTEND=noninteractive
|
||||||
|
apt-get update -q
|
||||||
|
apt-get install -y -q pigz git python3-venv curl unzip
|
||||||
|
snap install aws-cli --classic
|
||||||
|
S3="/snap/bin/aws s3"
|
||||||
|
( while sleep 300; do upload_log; done ) &
|
||||||
|
|
||||||
|
mkdir -p /work && cd /work
|
||||||
|
$S3 cp "s3://$BUCKET/$RUN/repo.tar.gz" . --only-show-errors
|
||||||
|
mkdir repo && tar -xzf repo.tar.gz -C repo && cd repo
|
||||||
|
curl -LsSf https://astral.sh/uv/install.sh | sh
|
||||||
|
export PATH="$HOME/.local/bin:$PATH"
|
||||||
|
uv venv -q --python 3.12 .venv && uv pip install -q -e ".[dev]"
|
||||||
|
mkdir -p refs && git clone -q --depth 1 https://github.com/Illumina/BeadArrayFiles.git refs/BeadArrayFiles
|
||||||
|
|
||||||
|
SRC=s3://1000genomes/release/20130502/supporting/hd_genotype_chip/broad_intensities
|
||||||
|
D=/work/omni && mkdir -p $D/data
|
||||||
|
for f in HumanOmni2.5-4v1-Multi_B.bpm HumanOmni2.5-4v1-Multi_B.egt HumanOmni2.5-4v1_B.csv; do
|
||||||
|
$S3 cp --no-sign-request "$SRC/$f" $D/ --only-show-errors
|
||||||
|
done
|
||||||
|
echo "=== streaming 107 GB IDAT+GTC archive (no copy of the .tgz is kept) $(date -u +%T) ==="
|
||||||
|
$S3 cp --no-sign-request "$SRC/Omni25_idats_gtcs_2141_samples.tgz" - | pigz -dc | tar -x -C $D/data
|
||||||
|
echo "=== extracted $(du -sh $D/data | cut -f1) at $(date -u +%T) ==="
|
||||||
|
|
||||||
|
echo "=== validation $(date -u +%T) ==="
|
||||||
|
.venv/bin/python -m pytest tests/unit -q -p no:cacheprovider
|
||||||
|
.venv/bin/python scripts/validate_p1_large.py --root $D/data --bpm $D/HumanOmni2.5-4v1-Multi_B.bpm \
|
||||||
|
--egt $D/HumanOmni2.5-4v1-Multi_B.egt --csv $D/HumanOmni2.5-4v1_B.csv --out /work/summary.json \
|
||||||
|
--workers "$(nproc)" || echo "validation reported differences (see summary.json)"
|
||||||
|
$S3 cp /work/summary.json "s3://$BUCKET/$RUN/summary.json" --only-show-errors
|
||||||
|
echo "=== done $(date -u +%T) ==="
|
||||||
@@ -0,0 +1,225 @@
|
|||||||
|
#!/usr/bin/env python3
|
||||||
|
"""Phase P1 validation on a large data set (e.g. 1000 Genomes Omni2.5, 2141 samples with vendor GTCs).
|
||||||
|
|
||||||
|
For every GTC: parse, compare every field with Illumina BeadArrayFiles (BSD-2) and test the byte-identical
|
||||||
|
round-trip. For every IDAT pair: parse, pair-check and compare the address set with the manifest. BPM/EGT/CSV are
|
||||||
|
compared with BeadArrayFiles once. Writes a JSON summary; exit code 0 only if everything matched.
|
||||||
|
|
||||||
|
Usage: validate_p1_large.py --root DIR --bpm F --egt F [--csv F] --out summary.json [--workers N] [--limit N]
|
||||||
|
"""
|
||||||
|
|
||||||
|
from __future__ import annotations
|
||||||
|
|
||||||
|
import argparse
|
||||||
|
import json
|
||||||
|
import os
|
||||||
|
import sys
|
||||||
|
import time
|
||||||
|
import traceback
|
||||||
|
from concurrent.futures import ProcessPoolExecutor
|
||||||
|
from pathlib import Path
|
||||||
|
from typing import Any
|
||||||
|
|
||||||
|
import numpy as np
|
||||||
|
|
||||||
|
ROOT = Path(__file__).resolve().parents[1]
|
||||||
|
sys.path.insert(0, str(ROOT / "refs" / "BeadArrayFiles"))
|
||||||
|
|
||||||
|
|
||||||
|
def check_gtc(path: str) -> dict[str, Any]:
|
||||||
|
from module.GenotypeCalls import GenotypeCalls # type: ignore[import-not-found]
|
||||||
|
|
||||||
|
from simurg.formats.gtc import encode_gtc, read_gtc
|
||||||
|
|
||||||
|
res: dict[str, Any] = {"file": path, "ok": True, "diffs": []}
|
||||||
|
try:
|
||||||
|
g = read_gtc(path, check_write_complete=False)
|
||||||
|
raw = Path(path).read_bytes()
|
||||||
|
res["version"] = g.version
|
||||||
|
res["autocall_version"] = g.autocall_version
|
||||||
|
res["roundtrip_identical"] = encode_gtc(g) == raw
|
||||||
|
if not res["roundtrip_identical"]:
|
||||||
|
res["roundtrip_canonical_identical"] = encode_gtc(g, canonical=True) == raw
|
||||||
|
res["size_ours_vs_file"] = [len(encode_gtc(g)), len(raw)]
|
||||||
|
b = GenotypeCalls(path, check_write_complete=False)
|
||||||
|
pairs = {
|
||||||
|
"raw_x": (g.raw_x, b.get_raw_x_intensities()),
|
||||||
|
"raw_y": (g.raw_y, b.get_raw_y_intensities()),
|
||||||
|
"genotypes": (g.genotypes, np.array(b.get_genotypes(), dtype=np.uint8)),
|
||||||
|
"scores": (g.genotype_scores, b.get_genotype_scores()),
|
||||||
|
}
|
||||||
|
if g.version >= 4:
|
||||||
|
pairs["baf"] = (g.b_allele_freqs, b.get_ballele_freqs())
|
||||||
|
pairs["lrr"] = (g.logr_ratios, b.get_logr_ratios())
|
||||||
|
for k, (a, c) in pairs.items():
|
||||||
|
if not np.array_equal(np.asarray(a), np.asarray(c), equal_nan=a.dtype.kind == "f"):
|
||||||
|
res["diffs"].append(k)
|
||||||
|
scalars = {
|
||||||
|
"call_rate": (g.call_rate, b.get_call_rate()),
|
||||||
|
"gender": (g.gender, b.get_gender()),
|
||||||
|
"logr_dev": (g.logr_dev, b.get_logr_dev()),
|
||||||
|
"gc10": (g.gc10, b.get_gc10()),
|
||||||
|
"gc50": (g.gc50, b.get_gc50()),
|
||||||
|
"num_calls": (g.num_calls, b.get_num_calls()),
|
||||||
|
"sample_name": (g.sample_name, b.get_sample_name()),
|
||||||
|
"base_calls": (g.base_call_strings(), [x.decode() for x in b.get_base_calls()]),
|
||||||
|
}
|
||||||
|
for k, (a, c) in scalars.items():
|
||||||
|
same = (a == c) or (isinstance(a, float) and np.isnan(a) and np.isnan(c))
|
||||||
|
if not same:
|
||||||
|
res["diffs"].append(k)
|
||||||
|
xf_b = [
|
||||||
|
(t.offset_x, t.offset_y, t.scale_x, t.scale_y, t.shear, t.theta) for t in b.get_normalization_transforms()
|
||||||
|
]
|
||||||
|
xf_g = [
|
||||||
|
tuple(float(x[k]) for k in ("offset_x", "offset_y", "scale_x", "scale_y", "shear", "theta"))
|
||||||
|
for x in g.normalization_transforms
|
||||||
|
]
|
||||||
|
if [tuple(float(v) for v in t) for t in xf_b] != xf_g:
|
||||||
|
res["diffs"].append("normalization_transforms")
|
||||||
|
res["ok"] = not res["diffs"]
|
||||||
|
except Exception as e: # noqa: BLE001 - report everything
|
||||||
|
res.update(ok=False, error=f"{type(e).__name__}: {e}", trace=traceback.format_exc(limit=3))
|
||||||
|
return res
|
||||||
|
|
||||||
|
|
||||||
|
def check_idat_pair(args: tuple[str, str]) -> dict[str, Any]:
|
||||||
|
from simurg.formats.idat import check_pair, read_idat
|
||||||
|
|
||||||
|
grn, red = args
|
||||||
|
res: dict[str, Any] = {"grn": grn, "ok": True}
|
||||||
|
try:
|
||||||
|
g, r = read_idat(grn), read_idat(red)
|
||||||
|
check_pair(g, r)
|
||||||
|
res.update(
|
||||||
|
n=g.n_addresses,
|
||||||
|
chip=g.chip_type,
|
||||||
|
midblock_itemsize=g.mid_block_itemsize,
|
||||||
|
scan=g.scan_time_iso,
|
||||||
|
name_ok=(g.barcode in Path(grn).name),
|
||||||
|
)
|
||||||
|
except Exception as e: # noqa: BLE001
|
||||||
|
res.update(ok=False, error=f"{type(e).__name__}: {e}")
|
||||||
|
return res
|
||||||
|
|
||||||
|
|
||||||
|
def main() -> int:
|
||||||
|
ap = argparse.ArgumentParser()
|
||||||
|
ap.add_argument("--root", required=True, type=Path)
|
||||||
|
ap.add_argument("--bpm", required=True, type=Path)
|
||||||
|
ap.add_argument("--egt", required=True, type=Path)
|
||||||
|
ap.add_argument("--csv", type=Path)
|
||||||
|
ap.add_argument("--out", required=True, type=Path)
|
||||||
|
ap.add_argument("--workers", type=int, default=os.cpu_count() or 1)
|
||||||
|
ap.add_argument("--limit", type=int, default=0)
|
||||||
|
a = ap.parse_args()
|
||||||
|
|
||||||
|
from simurg.formats.bpm import read_bpm
|
||||||
|
from simurg.formats.discovery import walk_files
|
||||||
|
from simurg.formats.egt import read_egt
|
||||||
|
|
||||||
|
t0 = time.time()
|
||||||
|
summary: dict[str, Any] = {"root": str(a.root), "started": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime())}
|
||||||
|
|
||||||
|
gtcs = [str(p) for p in walk_files(a.root, (".gtc",))]
|
||||||
|
idats = walk_files(a.root, (".idat", ".idat.gz"))
|
||||||
|
if a.limit:
|
||||||
|
gtcs = gtcs[: a.limit]
|
||||||
|
grn = {p.name.lower().replace("_grn.idat", ""): p for p in idats if "_grn.idat" in p.name.lower()}
|
||||||
|
red = {p.name.lower().replace("_red.idat", ""): p for p in idats if "_red.idat" in p.name.lower()}
|
||||||
|
pairs = [(str(grn[k]), str(red[k])) for k in sorted(set(grn) & set(red))]
|
||||||
|
if a.limit:
|
||||||
|
pairs = pairs[: a.limit]
|
||||||
|
summary["n_gtc"] = len(gtcs)
|
||||||
|
summary["n_idat_pairs"] = len(pairs)
|
||||||
|
summary["unpaired_idats"] = len(set(grn) ^ set(red))
|
||||||
|
|
||||||
|
with ProcessPoolExecutor(a.workers) as ex:
|
||||||
|
gres = list(ex.map(check_gtc, gtcs, chunksize=4))
|
||||||
|
ires = list(ex.map(check_idat_pair, pairs, chunksize=8))
|
||||||
|
summary["gtc_ok"] = sum(r["ok"] for r in gres)
|
||||||
|
summary["gtc_roundtrip_identical"] = sum(bool(r.get("roundtrip_identical")) for r in gres)
|
||||||
|
summary["gtc_versions"] = sorted({r.get("version") for r in gres if "version" in r})
|
||||||
|
summary["gtc_autocall_versions"] = sorted({r.get("autocall_version") for r in gres if "autocall_version" in r})
|
||||||
|
summary["gtc_failures"] = [r for r in gres if not r["ok"]][:50]
|
||||||
|
summary["gtc_roundtrip_failures"] = [r for r in gres if r.get("roundtrip_identical") is False][:20]
|
||||||
|
summary["idat_ok"] = sum(r["ok"] for r in ires)
|
||||||
|
summary["idat_failures"] = [r for r in ires if not r["ok"]][:50]
|
||||||
|
summary["idat_midblock_itemsizes"] = sorted({r.get("midblock_itemsize") for r in ires if r["ok"]})
|
||||||
|
summary["idat_chip_types"] = sorted({r.get("chip") for r in ires if r["ok"]})
|
||||||
|
|
||||||
|
# manifest / cluster file vs BeadArrayFiles
|
||||||
|
import module # type: ignore[import-not-found]
|
||||||
|
|
||||||
|
m = read_bpm(a.bpm)
|
||||||
|
b = module.BeadPoolManifest(str(a.bpm))
|
||||||
|
summary["bpm"] = {
|
||||||
|
"num_loci": m.num_loci,
|
||||||
|
"names": m.names == b.names,
|
||||||
|
"snps": m.snps == b.snps,
|
||||||
|
"chroms": m.chroms == b.chroms,
|
||||||
|
"map_infos": m.map_infos.tolist() == b.map_infos,
|
||||||
|
"address_a": m.address_a.tolist() == b.addresses,
|
||||||
|
"assay_types": m.assay_types.tolist() == b.assay_types,
|
||||||
|
"normalization_lookups": m.normalization_lookups().tolist() == b.normalization_lookups,
|
||||||
|
"n_bins": int(len(np.unique(m.norm_bin_ids()))),
|
||||||
|
}
|
||||||
|
e = read_egt(a.egt)
|
||||||
|
with open(a.egt, "rb") as fh:
|
||||||
|
be = module.ClusterFile.read_cluster_file(fh)
|
||||||
|
bad = 0
|
||||||
|
|
||||||
|
def same(a: float, b: float) -> bool:
|
||||||
|
return bool(np.float32(a) == np.float32(b) or (np.isnan(a) and np.isnan(b)))
|
||||||
|
|
||||||
|
for i, name in enumerate(e.names):
|
||||||
|
rec = be.get_record(name)
|
||||||
|
st = (rec.aa_cluster_stats, rec.ab_cluster_stats, rec.bb_cluster_stats)
|
||||||
|
if not all(
|
||||||
|
same(e.theta_mean[i, k], s.theta_mean)
|
||||||
|
and same(e.r_mean[i, k], s.r_mean)
|
||||||
|
and same(e.theta_dev[i, k], s.theta_dev)
|
||||||
|
and same(e.r_dev[i, k], s.r_dev)
|
||||||
|
and e.n[i, k] == s.N
|
||||||
|
for k, s in enumerate(st)
|
||||||
|
):
|
||||||
|
bad += 1
|
||||||
|
summary["egt"] = {
|
||||||
|
"num_records": e.num_records,
|
||||||
|
"dbv": e.data_block_version,
|
||||||
|
"extra_floats": e.extra_per_record,
|
||||||
|
"record_mismatches": bad,
|
||||||
|
}
|
||||||
|
if a.csv:
|
||||||
|
from simurg.formats.csv_manifest import compare_manifests, read_csv_manifest
|
||||||
|
|
||||||
|
summary["bpm_csv_problems"] = compare_manifests(m, read_csv_manifest(a.csv, load_sequences=False))
|
||||||
|
if gtcs:
|
||||||
|
g0 = __import__("simurg.formats.gtc", fromlist=["read_gtc"]).read_gtc(gtcs[0], check_write_complete=False)
|
||||||
|
summary["gtc_vs_bpm_num_snps"] = [g0.num_snps, m.num_loci]
|
||||||
|
summary["gtc_transforms_vs_bins"] = [len(g0.normalization_transforms), summary["bpm"]["n_bins"]]
|
||||||
|
|
||||||
|
summary["seconds"] = round(time.time() - t0, 1)
|
||||||
|
passed = (
|
||||||
|
summary["gtc_ok"] == len(gtcs)
|
||||||
|
and summary["idat_ok"] == len(pairs)
|
||||||
|
and all(v for k, v in summary["bpm"].items() if isinstance(v, bool))
|
||||||
|
and summary["egt"]["record_mismatches"] == 0
|
||||||
|
and not summary.get("bpm_csv_problems")
|
||||||
|
)
|
||||||
|
summary["passed"] = passed
|
||||||
|
a.out.parent.mkdir(parents=True, exist_ok=True)
|
||||||
|
a.out.write_text(json.dumps(summary, indent=2, default=str) + "\n")
|
||||||
|
print(
|
||||||
|
json.dumps(
|
||||||
|
{
|
||||||
|
k: summary[k]
|
||||||
|
for k in ("n_gtc", "gtc_ok", "gtc_roundtrip_identical", "n_idat_pairs", "idat_ok", "passed", "seconds")
|
||||||
|
}
|
||||||
|
)
|
||||||
|
)
|
||||||
|
return 0 if passed else 1
|
||||||
|
|
||||||
|
|
||||||
|
if __name__ == "__main__":
|
||||||
|
sys.exit(main())
|
||||||
Reference in new issue
Block a user