← DESeq Desk / API
Tokens

Drive DESeq Desk from your own code

A review of the statistics, not of the biology. The model reads what your browser (or your script) computed from the counts; it never sees the counts, never states what a gene does and never recomputes a number. "Sound" means the design is replicated and nothing flagged at medium or high undermines the calls - not that a gene has any function, and a gene without an adjusted p-value is not evidence of no change.

Everything the web page does is available over HTTP. Run the analysis with the page's own deseq.js and deseqkit.js (pydeseq2 0.5.4 conventions: the 10-count gene filter, median-of-ratios size factors, gene-wise, trend and MAP dispersions, IRLS fold changes, Cook's outlier replacement and refit, the Wald test, Cook's and independent filtering with Benjamini-Hochberg adjustment and optional LFC shrinkage, minimised with a port of scipy's L-BFGS-B), send the facts, and get back a verdict (sound, caveated, unreliable) and either a review of every metric and of the result or a pydeseq2 script that reproduces it and runs follow-up checks. The natural loop: load the counts, review, script, change the design or a filter, re-check. DESeq Desk is derived from the agent skill @k-dense-ai/pydeseq2 (k-dense-ai/scientific-agent-skills, skill author K-Dense Inc.) and its scripts/run_deseq2_analysis.py.

Two lanes: the task field

taskreturnsextra input
reviewA reading of every metric, what the result shows (how many genes, which direction, the strongest calls by id), changes to try (min_counts, a batch term, leaving out a named sample, shrinkage, a thresholded test, alpha, a VST PCA, more replicates), your claims judged against the facts, a methods paragraph stating the exact design, contrast, filter and alpha, and what the result cannot show.none
scriptThe fixes and one complete Python script: COUNTS and METADATA read with pd.read_csv(..., index_col=0), the same alignment, reference-first categories and count filter, DeseqDataSet and DeseqStats with every browser setting as a literal, lfc_shrink exactly when the browser shrank, an EXPECTED dict of browser values checked with math.isclose, then the fixes, all inside main().decision: the text of an earlier review run (optional)

Input fields

fieldrequiredmeaning
taskyes"review" or "script". An unknown value makes the model pick the closest lane and say which in lane.
factsyesA JSON-encoded string holding the browser's DESeq2 run - see below. The page builds it with DeseqKit.buildInput. The counts themselves are never sent.
titlenoA label for the analysis.
contextnoYour notes: the experiment, the samples, what you want to conclude (the page sends up to 2,500 characters).
questionnoAnswered in the first summary bullet (up to 1,500 characters).
decisionnoScript lane only: the text of an earlier review (up to 5,000 characters).
retry_notenoSent by the page on its one reformat retry after a malformed reply.

The facts string

facts is a JSON string, not an object. It holds settings (pydeseq2_version 0.5.4; design, e.g. ~batch + condition; condition_column, reference, test and contrast; batch_column; levels with the reference first; design_columns; min_counts; alpha; refit_cooks, cooks_filter, independent_filter; lfc_shrink and shrink_coeff; the file names, separators, orientation and any skipped featureCounts columns; runner_equivalent and runner_differences; for synthetic data synthetic), samples (sample, condition, batch, size factor, library size), metrics (M1..: samples, genes tested, size-factor range, significant genes up and down, significant with |log2FC| > 1, the dispersion trend, the prior variance, dispersion outliers, Cook's outliers, the independent-filtering threshold, the p-value histogram, the shrinkage prior and, for synthetic data, true and false positives), flags (F1.., severity high / medium / low), top_genes (up to 20 by padj), significant_total, browser_verdict and expected - the values a pydeseq2 reproduction must match: sample and gene counts, the significant counts, every size factor, the trend coefficients, the prior variance and baseMean, log2FoldChange, pvalue and padj for the top five genes.

To copy the facts without writing code, run DESeq2 on the page and press Report .json: the file carries the exact facts object under facts, next to the full results arrays. A finished result's Download .json carries it too, under browser. Send it back as a string: json.dumps(facts), JSON.stringify(facts) or your language's equivalent.

Building the body

To build a body that matches the page byte for byte, run the page's own modules in Node 18 or later: save detmath.js, lbfgsb.js, deseq.js, ioparse.js and deseqkit.js next to the script below. All five export themselves with module.exports. 

// make-body.js - build the exact body the page sends, with the page's own code.
//   node make-body.js counts.csv samples.csv condition control treated review "Title" "notes" > body.json
const fs = require("fs");
const IO = require("./ioparse.js"), D = require("./deseq.js"), K = require("./deseqkit.js");
const [countsFile, sheetFile, condCol, ref, test, lane = "review", title = "", context = "", batchCol = ""] = process.argv.slice(2);
const sheet = IO.parseSheet(fs.readFileSync(sheetFile, "utf8"));
const counts = IO.parseCounts(fs.readFileSync(countsFile, "utf8"), sheet.ids);
const al = IO.align(counts, sheet, condCol, batchCol || null);
const s = { condCol, batchCol, ref, test, minCounts: 10, alpha: 0.05, shrink: true, cooksFilter: true, independentFilter: true, refitCooks: true };
const res = D.run(Object.assign({ genes: al.genes, counts: al.counts, samples: al.samples, cond: al.cond, batch: al.batch }, s, { batchCol: batchCol || "batch" }));
const X = K.analyze(res, { al, settings: s, title, duplicateGenes: IO.duplicates(al.genes),
  parsed: { countsName: countsFile, sheetName: sheetFile, sampleIdColumn: sheet.idColumn, orientation: counts.orientation,
            featureCounts: counts.featureCounts, skipped: counts.skipped, countsSep: counts.delim, sheetSep: sheet.delim } });
const body = K.mustBeObject(K.buildInput(X, { lane, title, context }));
fs.writeFileSync("deseq2_results.csv", K.resultsCsv(X, false));   // pydeseq2 results_df layout
console.error("browser verdict:", X.hint, "| tested:", X.n, "| significant:", X.nSig, "| flags:", X.flags.map(f => f.id + " " + f.category).join(", "));
console.error("idempotency key: deseq-desk:" + body.task + ":" + K.hashInput(body) + ":a1");
process.stdout.write(JSON.stringify(body));
# Or build the body in any language from a facts object you already hold, for example the
# "facts" key of the page's "Report .json" download. facts must go in as a STRING.
import json

report = json.load(open("deseq-desk-knockdown-vs-control-batch-adjusted-report.json"))   # the page's download
body = {
    "task": "review",
    "title": "Knockdown vs control, batch-adjusted",
    "context": "Four knockdown and four control cultures, sequenced in two batches ...",
    "question": "Is adjusting for batch the right design here?",
    "facts": json.dumps(report["facts"], separators=(",", ":")),
}
json.dump(body, open("body.json", "w"))

Base URL and the envelope

Every endpoint lives under https://api.skillsafe.ai/v1/app-api and every response uses the same envelope, so one helper covers the whole API: 

{"ok": true, "data": {"job_id": "job_...", "status": "queued"}}
{"ok": false, "error": {"code": "payment_required", "message": "..."}}

The token is minted for this app (the guest endpoint takes {"slug":"deseq-desk"} in its body), so no slug header is needed afterwards. Send it as Authorization: Bearer ....

The input object IS the request body. There is no {"input": ...} wrapper. A wrapped body is answered with an unknown field 'input' warning, and the model never sees your text.

Error codes

statuscodewhat to do
400validation_errorA field is missing or the wrong type. Every field is a string: facts must be a JSON-encoded string, not an object.
401unauthorizedThe token is missing, malformed or expired. Get a new one from the token page.
402payment_requiredThe balance is below min_credits. Call /estimate first and top up.
403forbiddenThe token is valid but not for this app, or a guest token tried a metered run. A guest cannot run; sign in for a personal token.
404not_foundUnknown job id, or the app slug does not exist.
409conflictThe same Idempotency-Key was replayed with a different body. Change the key or send the original input.
429rate_limitedToo many requests. Back off and retry; do not tight-loop.
5xxinternalA server-side failure. Retry with the SAME Idempotency-Key so you are not billed twice.

1. A tiny client

One helper that sends the token, unwraps data and raises on ok: false. The token comes from the token page (Copy token or Copy shell export); step 2 covers the kinds of token and minting one from code. 

# Every call is the same three things: the base URL, your bearer token,
# and a JSON body. Keep the token in a shell variable.
BASE="https://api.skillsafe.ai/v1/app-api"
SLUG="deseq-desk"
TOKEN="$SKILLSAFE_TOKEN"   # from https://deseq-desk.skillsafe.ai/tokens.html

call() {                  # call <path> [json-body]
  if [ -n "$2" ]; then
    curl -sS -X POST "$BASE/$1" \
      -H "Authorization: Bearer $TOKEN" \
      -H "Content-Type: application/json" \
      -d "$2"
  else
    curl -sS "$BASE/$1" -H "Authorization: Bearer $TOKEN"
  fi
}

2. Get a token

The easiest route is the token page: it shows the token this browser already holds, with Copy token and Copy shell export buttons, and a sign-in button for a personal token. A guest token, minted with POST /guest and {"slug":"deseq-desk"}, can call /me and /estimate; the run is metered, so /run and /run-stream need a personal token. 

# The token page is the shortest path. It shows the token this browser holds and
# hands you a ready-made shell export:
#
#   https://deseq-desk.skillsafe.ai/tokens.html
#   export SKILLSAFE_TOKEN="..."
#
# To mint a guest token from the command line instead. A guest token is enough
# for /me and /estimate; a run needs a personal token from signing in.
curl -sS -X POST "https://api.skillsafe.ai/v1/app-api/guest" \
  -H "Content-Type: application/json" -d '{"slug":"deseq-desk"}'
# {"ok":true,"data":{"token":"...","subject_type":"guest"}}

3. Check the session and the balance

call me
# {"ok":true,"data":{"subject_type":"user","username":"you","credits":51234}}

4. Price the run (free)

/estimate returns the model binding and the credits a run would reserve. It creates no job and charges nothing. Expect model_alias gpt-terra and markup_bps 1000 (a 10% markup). hold_credits is a reservation, not the price: it is held against your balance while the run executes and released afterwards. min_credits is the least balance that can start a run. What you actually pay is charged_credits, reported on the finished job and in the done event, and it is usually far lower than the hold. The body is the input object itself, with no {"input": ...} wrapper. /estimate does not validate the body, so check the shape yourself: an object whose every value is a string, task equal to review or script, facts non-empty, and facts a JSON string that parses to an object (this is what the page's own guard, DeseqKit.mustBeObject, refuses to spend without). 

# body.json is the input object itself - no {"input": ...} wrapper. Build it with
# make-body.js above, or by hand. estimate does not validate it, so check the shape first:
python3 -c 'import json;b=json.load(open("body.json"));assert isinstance(b,dict) and b.get("task") in ("review","script") and all(isinstance(v,str) for v in b.values()) and all(b.get(k,"").strip() for k in ("facts",)) and isinstance(json.loads(b["facts"]),dict)'
INPUT=$(cat body.json)

call estimate "$INPUT"
# {"ok":true,"data":{"model":"...","model_alias":"gpt-terra",
#   "markup_bps":1000,"hold_credits":...,"min_credits":...,"sponsor_enabled":false,
#   "warnings":[]}}
#
# estimate creates no job and charges nothing. hold_credits is RESERVED, not the
# price; charged_credits after the run is the actual cost, usually far lower.

5. Run it, then poll

POST /run returns a job_id; poll GET /jobs/{id} until it is terminal. The reply is a string at data.output.output: JSON.parse it (step 7). Send an Idempotency-Key built from the lane, a hash of the input and the attempt number, deseq-desk:<lane>:<hash>:a<attempt> (for example deseq-desk:review:17z4spz1v70svb:a1), so a retried request returns the same job instead of billing a second run. Use one key per distinct input: changed counts, settings or notes (so changed facts) or a changed review are a new hash, the same data in the other lane is a new key, and replaying an old key with a different body is a 409. The page uses DeseqKit.hashInput(body) for the hash (it covers task, title, context, facts, decision and question; make-body.js prints the key); any stable digest of the body works from other languages. Leave retry_note out of the hash and bump the attempt instead. 

# Always send an Idempotency-Key derived from the input. A retried request with
# the same key returns the SAME job instead of billing a second run.
LANE=$(printf '%s' "$INPUT" | python3 -c 'import sys,json;print(json.load(sys.stdin)["task"])')   # review or script
KEY="deseq-desk:$LANE:$(printf '%s' "$INPUT" | shasum -a 256 | cut -c1-16):a1"

JOB=$(curl -sS -X POST "$BASE/run" \
  -H "Authorization: Bearer $TOKEN" \
  -H "Content-Type: application/json" \
  -H "Idempotency-Key: $KEY" \
  -d "$INPUT" | python3 -c 'import sys,json;print(json.load(sys.stdin)["data"]["job_id"])')

while :; do
  OUT=$(call "jobs/$JOB")
  STATUS=$(printf '%s' "$OUT" | python3 -c 'import sys,json;print(json.load(sys.stdin)["data"]["status"])')
  [ "$STATUS" = "succeeded" ] && break
  [ "$STATUS" = "failed" ] && echo "$OUT" && exit 1
  sleep 2
done

# {"ok":true,"data":{"job_id":"job_...","status":"succeeded",
#   "output":{"output":"{\"lane\":\"review\",\"verdict\":\"caveated\",\"headline\":\"...\", ...}"},
#   "charged_credits":...,"truncated":false}}
printf '%s' "$OUT" | python3 -c 'import sys,json;print(json.load(sys.stdin)["data"]["output"]["output"])' > reply.json

6. Or stream it

POST /run-stream takes the same body and headers and answers with server-sent events: job (the job id), delta (chunks of the reply) and done (the status, charged_credits, truncated and, when present, the full output). A browser page may receive only tick heartbeats and then done, never a delta, so take the reply from done.output.output when it is there, fall back to the concatenated deltas, and fall back again to GET /jobs/{id}. 

# Server-sent events. `delta` events carry chunks of the reply; `done` carries the
# status, charged_credits and the truncated flag. Ignore `tick` heartbeats.
curl -N -X POST "$BASE/run-stream" \
  -H "Authorization: Bearer $TOKEN" \
  -H "Content-Type: application/json" \
  -H "Idempotency-Key: $KEY" \
  -H "Accept: text/event-stream" \
  -d "$INPUT"

# event: job    {"job_id":"job_..."}
# event: delta  {"text":"{\"lane\":\"review\",\"verdict\":\"caveated\",\"headline\":\"The"}
# event: done   {"status":"succeeded","charged_credits":...,"truncated":false}

7. Parse the reply

The reply is a JSON object serialised as a string. Parse it, then check the lane.

# The reply is a JSON string inside data.output.output. Pull it out and parse it:
printf '%s' "$JOB" | python3 -c 'import sys,json;r=json.loads(json.load(sys.stdin)["output"]["output"]);print(r["verdict"],r["headline"])'

Invariants worth asserting

Worked example: review

A real request: the synthetic knockdown experiment from the page's first example (4 knockdown and 4 control samples in two batches, design ~batch + condition, shrunk fold changes). The browser tests 2,854 genes, calls 133 significant and raises one low flag, so its read is sound. The body, with facts abbreviated (send the full string from make-body.js or the page's report): 

{
 "task": "review",
 "title": "Knockdown vs control, batch-adjusted",
 "context": "Four knockdown and four control cultures, sequenced in two batches with both conditions in each batch. I want to say the knockdown changes a few hundred genes and report the strongest ones.",
 "facts": "{\"settings\":{\"pydeseq2_version\":\"0.5.4\",\"design\":\"~batch + condition\",\"contrast\":[\"condition\",\"knockdown\",\"control\"],\"levels\":[\"control\",\"knockdown\"],\"min_counts\":10,\"alpha\":0.05,\"lfc_shrink\":true,\"shrink_coeff\":\"condition[T.knockdown]\",\"counts_file\":\"counts.csv\",\"sample_sheet_file\":\"samples.csv\"},\"samples\":[{\"sample\":\"control_1\",\"condition\":\"control\",\"size_factor\":0.886628,\"library_size\":925418,\"batch\":\"b1\"},{\"sample\":\"control_2\",\"condition\":\"control\",\"size_factor\":1.26242,\"library_size\":1349309,\"batch\":\"b2\"},\"... 6 more\"],\"metrics\":[{\"id\":\"M1\",\"metric\":\"samples\",\"value\":[4,4],\"basis\":\"samples in the contrast: control (reference) then knockdown; all levels: control 4, knockdown 4; 8 samples, 3 design coefficients, 5 residual degrees of freedom\",\"residual_df\":5},{\"id\":\"M2\",\"metric\":\"genes_tested\",\"value\":2854,\"basis\":\"genes with at least 10 counts in total, of 3000 in the file (the script's --min-counts filter)\",\"fraction\":0.9513333333333334},\"... M3 to M13\"],\"flags\":[{\"id\":\"F1\",\"severity\":\"low\",\"category\":\"dispersion\",\"message\":\"The dispersion prior variance sits at its 0.25 floor.\",\"refs\":\"M7\"}],\"top_genes\":[{\"gene\":\"GENE1253\",\"baseMean\":16186.9,\"log2FoldChange\":3.9054,\"lfcSE\":0.1535,\"stat\":25.532,\"pvalue\":8.76e-144,\"padj\":2.306e-140,\"log2FoldChange_mle\":3.9167},\"... 19 more\"],\"significant_total\":133,\"browser_verdict\":\"sound\",\"expected\":{\"n_samples\":8,\"n_genes_tested\":2854,\"n_significant\":133,\"n_up\":73,\"...\":\"32 more\"},\"expected_count\":36}",
 "question": "Is adjusting for batch the right design here?"
}

The reply the stand-in model gave during validation, which the page replays for free and which reconciles with zero disagreements:

{
  "lane": "review",
  "verdict": "sound",
  "headline": "On this synthetic dataset, the batch-adjusted design calls 133 of 2854 tested genes significant (73 up, 60 down), recovering 128 of 297 known true changes with 5 false positives.",
  "tldr": [
    "Answer: Yes - both batches (b1, b2) contain both control and knockdown samples, so condition is not confounded with batch, and the ~batch + condition design correctly adjusts for it.",
    "133 of the 2854 tested genes are significant at padj < 0.05 (73 up, 60 down in knockdown vs control) - a real but smaller call set than 'a few hundred'.",
    "This run is on synthetic data (seed 20260927, 313 true DE genes among 3000 genes in the file), so these numbers describe pipeline recovery, not a real knockdown experiment.",
    "The only browser flag (F1, low severity) notes the dispersion prior variance sits at its 0.25 floor; nothing is flagged medium or high, so the fit reads as sound.",
    "The strongest calls by padj are GENE1253 (log2FoldChange 3.91, padj 2.31e-140), GENE1667 (log2FoldChange 2.98, padj 2.07e-60) and GENE0757 (log2FoldChange 2.43, padj 1.49e-57), all higher in knockdown."
  ],
  "metrics": [
    {
      "id": "M1",
      "reading": "Four control and four knockdown samples (8 total) were fit with 3 design coefficients (Intercept, batch[T.b2], condition[T.knockdown]), leaving 5 residual degrees of freedom."
    },
    {
      "id": "M2",
      "reading": "2854 of the 3000 genes in the file passed the min_counts=10 filter (about 95.1%) and were tested."
    },
    {
      "id": "M3",
      "reading": "Median-of-ratios size factors across the 8 samples range from 0.733 to 1.26, a max/min ratio of about 1.72 - a moderate, not extreme, depth spread."
    },
    {
      "id": "M4",
      "reading": "133 genes are called significant at padj < 0.05 (about 4.66% of the 2854 tested genes): 73 up and 60 down in knockdown vs control."
    },
    {
      "id": "M5",
      "reading": "Of the 133 significant genes, 92 also have a shrunk |log2FoldChange| greater than 1."
    },
    {
      "id": "M6",
      "reading": "The dispersion trend is parametric (a0=0.0253, a1=2.34), fit over 2241 genes in 3 iterations - an ordinary fit type rather than a mean fallback."
    },
    {
      "id": "M7",
      "reading": "The dispersion prior variance is 0.25, at pydeseq2's floor (the raw estimate - squared MAD of log residuals 0.5518 minus trigamma(5/2) - would be lower); this is what F1 flags."
    },
    {
      "id": "M8",
      "reading": "22 genes (about 0.771% of tested genes) have gene-wise dispersion more than 2 prior SDs above the trend and keep their gene-wise value instead of the MAP value."
    },
    {
      "id": "M9",
      "reading": "No genes were flagged as Cook's outliers (threshold F(0.99;3,5)=12.06), and since no design cell has 7 or more replicates, none were refit."
    },
    {
      "id": "M10",
      "reading": "Independent filtering set a baseMean threshold of 4.01 (the 7.76th percentile), removing 222 genes (about 7.78% of tested genes) that had a p-value from getting a padj."
    },
    {
      "id": "M11",
      "reading": "The p-value histogram across ten 0.1-wide bins (391, 251, 255, 280, 278, 264, 266, 311, 267, 291) shows a spike near zero and is otherwise roughly flat, consistent with genuine signal rather than a hill shape; 284 genes had a raw p-value below 0.05 before adjustment."
    },
    {
      "id": "M12",
      "reading": "LFC shrinkage used a Cauchy prior scale of 0.283 for condition[T.knockdown]; the reported log2FoldChange and lfcSE are shrunk while stat and pvalue come from the unshrunk fit."
    },
    {
      "id": "M13",
      "reading": "On the synthetic truth set, 128 of the 297 truly changed genes tested were recovered at padj < 0.05, with 5 false positives (about 3.76% false discovery share) and 169 true changes missed."
    }
  ],
  "findings": [
    "This run is on synthetic data (seed 20260927; 313 of 3000 genes truly changed), so these counts describe pipeline recovery rather than a biological experiment.",
    "133 of 2854 tested genes are significant at padj < 0.05 (73 up, 60 down in knockdown vs control), about 4.66% of tested genes; 92 of those also have |log2FoldChange| > 1 (shrunk).",
    "The strongest calls by padj are GENE1253 (log2FoldChange 3.91, padj 2.31e-140), GENE1667 (log2FoldChange 2.98, padj 2.07e-60) and GENE0757 (log2FoldChange 2.43, padj 1.49e-57), all higher in knockdown.",
    "Size factors range from 0.733 (knockdown_2) to 1.26 (control_2), a max/min of about 1.72, and library sizes (863,215 to 1,349,309) show no sample an order of magnitude off the rest.",
    "Against the 297 truly changed genes tested (of 313 total), the pipeline recovered 128 with 5 false positives and missed 169."
  ],
  "parameters": [
    {
      "change": "Run a variance-stabilized PCA (dds.vst(), then a PCA of dds.layers['vst_counts']) across the 8 samples.",
      "why": "A suggestion to test whether the two batches hide any structure beyond what the ~batch + condition design already models; not required by any flag here.",
      "refs": ""
    },
    {
      "change": "Re-run DeseqStats with lfc_null=1 and alt_hypothesis='greaterAbs' instead of filtering the shrunk log2FoldChange post hoc.",
      "why": "A suggestion to test: this directly asks which genes clear a fold-change floor with a formal test, closer to a 'genes that changed' claim than filtering the 133 significant calls by |log2FoldChange| > 1 (92 pass) after the fact.",
      "refs": ""
    }
  ],
  "claims": [
    {
      "claim": "The knockdown changes a few hundred genes",
      "support": "partly",
      "why": "133 of the 2854 tested genes are significant at padj < 0.05 (73 up, 60 down) - a genuine, sizeable call set, but smaller than 'a few hundred' typically implies, and this run is on synthetic data rather than the described cultures."
    }
  ],
  "methods": "Differential expression was assessed with pydeseq2 0.5.4 using the design formula ~batch + condition on 8 samples (4 control, 4 knockdown), testing the contrast condition knockdown versus control (reference). Genes with fewer than 10 total counts were removed before fitting, leaving 2854 tested genes; size factors were estimated by the median-of-ratios method. Cook's distances were used for filtering and, in design cells with 7 or more replicates, for outlier refitting (refit_cooks=True), and independent filtering with Benjamini-Hochberg adjustment was applied at alpha=0.05 to call significance. Log2 fold changes and their standard errors were then shrunk with ds.lfc_shrink for the condition[T.knockdown] coefficient; reported test statistics and p-values are from the unshrunk fit.",
  "cautions": [
    "This run used synthetic data (seed 20260927; 313 of 3000 genes truly changed in the file), not a real culture experiment; treat these numbers as a pipeline check, not a biological result.",
    "padj < 0.05 controls the expected share of false calls in the significant set as a whole, not the chance that any one gene is a false call, and says nothing about what a called gene does.",
    "No gene lost a p-value to Cook's filtering here (0 outliers), but independent filtering removed padj for 222 genes with low baseMean (about 7.78% of tested genes); that is not evidence those genes are unchanged.",
    "The significant-gene list and the 92-gene |log2FoldChange| > 1 subset are statistical calls only; no gene function, pathway or cause is stated here - an enrichment analysis is the next step for biological interpretation."
  ],
  "next_steps": [
    "Confirm with the app whether the counts and sample sheet were meant to be the real culture experiment - settings.synthetic shows this run's data is synthetic (seed 20260927).",
    "Send the 133 significant genes (or the 92 with |log2FoldChange| > 1) to a pathway or enrichment tool for biological interpretation.",
    "Keep the ~batch + condition design - each batch (b1, b2) contains both control and knockdown samples, so batch is not confounded with condition.",
    "If a stricter effect-size cutoff is needed, rerun DeseqStats with lfc_null=1 and alt_hypothesis='greaterAbs' rather than filtering shrunk log2FoldChange values after the fact.",
    "Only chase the 222 genes removed by independent filtering (M10) if there is prior interest in those specific low-count genes; otherwise treat their exclusion as expected at this depth."
  ],
  "prescan_responses": [
    {
      "ref": "F1",
      "verdict": "confirmed",
      "note": "The dispersion prior variance is 0.25 and is at pydeseq2's floor (M7). This is a low-severity note about how much gene-wise dispersions can be pulled toward the trend; the dispersion trend, outlier counts and p-value histogram otherwise look ordinary, so it does not change the verdict."
    }
  ]
}

Worked example: script

The same experiment handed to the script lane with the review above as decision. The returned script was run with pydeseq2 0.5.4 on the page's counts.csv and samples.csv and printed no EXPECTED mismatch. Body (abbreviated as above, plus decision): 

{
 "task": "script",
 "title": "Knockdown vs control, batch-adjusted",
 "context": "I want a script that reproduces this analysis and then checks the samples for outliers with a VST PCA.",
 "facts": "{\"settings\":{\"pydeseq2_version\":\"0.5.4\",\"design\":\"~batch + condition\",\"contrast\":[\"condition\",\"knockdown\",\"control\"],\"levels\":[\"control\",\"knockdown\"],\"min_counts\":10,\"alpha\":0.05,\"lfc_shrink\":true,\"shrink_coeff\":\"condition[T.knockdown]\",\"counts_file\":\"counts.csv\",\"sample_sheet_file\":\"samples.csv\"},\"samples\":[{\"sample\":\"control_1\",\"condition\":\"control\",\"size_factor\":0.886628,\"library_size\":925418,\"batch\":\"b1\"},{\"sample\":\"control_2\",\"condition\":\"control\",\"size_factor\":1.26242,\"library_size\":1349309,\"batch\":\"b2\"},\"... 6 more\"],\"metrics\":[{\"id\":\"M1\",\"metric\":\"samples\",\"value\":[4,4],\"basis\":\"samples in the contrast: control (reference) then knockdown; all levels: control 4, knockdown 4; 8 samples, 3 design coefficients, 5 residual degrees of freedom\",\"residual_df\":5},{\"id\":\"M2\",\"metric\":\"genes_tested\",\"value\":2854,\"basis\":\"genes with at least 10 counts in total, of 3000 in the file (the script's --min-counts filter)\",\"fraction\":0.9513333333333334},\"... M3 to M13\"],\"flags\":[{\"id\":\"F1\",\"severity\":\"low\",\"category\":\"dispersion\",\"message\":\"The dispersion prior variance sits at its 0.25 floor.\",\"refs\":\"M7\"}],\"top_genes\":[{\"gene\":\"GENE1253\",\"baseMean\":16186.9,\"log2FoldChange\":3.9054,\"lfcSE\":0.1535,\"stat\":25.532,\"pvalue\":8.76e-144,\"padj\":2.306e-140,\"log2FoldChange_mle\":3.9167},\"... 19 more\"],\"significant_total\":133,\"browser_verdict\":\"sound\",\"expected\":{\"n_samples\":8,\"n_genes_tested\":2854,\"n_significant\":133,\"n_up\":73,\"...\":\"32 more\"},\"expected_count\":36}",
 "decision": "Verdict: sound.\nBatch-adjusted DESeq2 comparison of knockdown versus control calls 133 genes significant (73 up, 60 down) out of 2854 tested, a sound result with one low-severity flag.\n- 133 of 2854 tested genes are called significant at padj < 0.05, 73 up and 60 down in knockdown versus control (M4), about 4.66% of tested genes.\n- The strongest up call is GENE1253 (log2FoldChange 3.91, padj 2.31e-140, baseMean 16200), and the strongest down call by padj rank is GENE2444 (log2FoldChange -2.31, padj 2.33e-32).\n- Size factors range from 0.733 (knockdown_2) to 1.26 (control_2), a max/min of 1.72, an ordinary spread across the 8 samples (M3).\n- 92 of the 133 significant genes also pass |log2FoldChange| > 1 (M5), so most of the calls are large-effect.\n- On the synthetic truth set, 128 of 297 truly changed genes were recovered as significant with 5 false positives (M13), a 3.76% false discovery share among calls.\n- Try: Compare against a design without the batch term (~condition only) and see whether the significant gene count changes materially. (Directly tests whether adjusting for batch is doing real work here, answering the question of whether the batch term is warranted.)\n- Try: Run a thresholded test with lfc_null=1 and alt_hypothesis=\"greaterAbs\" in DeseqStats. (Focuses the call on the 92 genes with |log2FoldChange| > 1 (M5) rather than any nonzero change.)\n- Try: Run a variance-stabilized PCA (dds.vst(), then PCA of vst_counts) across the 8 samples. (Checks visually that samples separate by condition rather than by batch, supporting the ~batch + condition design.)\nNext steps:\n- Run a variance-stabilized PCA (dds.vst() then PCA) to confirm samples separate by condition rather than batch, supporting the batch + condition design.\n- Compare against a design without the batch term to see whether it changes the significant gene count.\n- Run a thresholded test with lfc_null=1 and alt_hypothesis=\"greaterAbs\" to focus on the 92 genes with |log2FoldChange| > 1 (M5).\n- Send the significant gene list to an enrichment or pathway analysis, since this result is statistical only.\n- Consider more replicates per condition in a follow-up experiment to raise the residual degrees of freedom above the current 5 (M1)."
}

Reply:

{
  "lane": "script",
  "verdict": "sound",
  "headline": "This script reproduces the batch-adjusted knockdown-vs-control pydeseq2 0.5.4 run (133 of 2854 genes significant, M4) and adds a VST PCA to check the samples for outliers or batch structure, as requested.",
  "tldr": [
    "The script reads counts.csv and samples.csv, re-fits the ~batch + condition pydeseq2 0.5.4 model exactly as the browser did (design_columns Intercept, batch[T.b2], condition[T.knockdown]), and checks the 36 EXPECTED values (n_significant 133, M4) with math.isclose tolerances.",
    "It also runs a VST PCA on the 8 samples, as you asked, to look for condition versus batch separation (4 control, 4 knockdown across batches b1/b2).",
    "Two more re-runs are included: dropping the batch term, and a thresholded test at lfc_null=1 focused on the 92 large-effect calls (M5).",
    "Verdict carried over from the review: sound, with one low-severity flag (F1) about the dispersion prior sitting at its 0.25 floor (M7)."
  ],
  "fixes": [
    {
      "fix": "Run a variance-stabilized PCA: dds.vst() then a PCA (via numpy SVD on the centered vst_counts) of the 8 samples, printing PC1/PC2 by sample.",
      "why": "Checks whether the 8 samples (4 control, 4 knockdown) separate by condition or by batch (b1, b2) on a variance-stabilized PCA, beyond the size factor spread already seen (M3: 0.733-1.26).",
      "refs": ""
    },
    {
      "fix": "Save the full results table to deseq2_results.csv and the padj < 0.05 subset to significant_genes.csv.",
      "why": "Keeps the full results table (2854 tested genes, M2) and just the 133 significant genes (M4) available for downstream use.",
      "refs": ""
    },
    {
      "fix": "Print the top 20 genes by padj (baseMean, log2FoldChange, pvalue, padj) to the console.",
      "why": "Gives a quick check against the review's top-gene table, for example GENE1253 (log2FoldChange 3.91, padj 2.31e-140).",
      "refs": ""
    },
    {
      "fix": "Save a volcano plot (log2FoldChange vs -log10 padj) as volcano.png using matplotlib's Agg backend.",
      "why": "Gives a visual read on the 2854 tested genes (M2) and where the 133 significant calls (M4) fall.",
      "refs": ""
    },
    {
      "fix": "Re-run DeseqStats with design=\"~condition\" (no batch term) and report how the significant count changes.",
      "why": "Tests whether the batch term (2 levels, b1/b2) is doing real work, against the 133 significant genes (M4) found with ~batch + condition.",
      "refs": ""
    },
    {
      "fix": "Re-run with a thresholded test: lfc_null=1.0 and alt_hypothesis=\"greaterAbs\" in DeseqStats, keeping the same design, and report the new significant count.",
      "why": "Focuses the call on the 92 genes already passing |log2FoldChange| > 1 (M5) rather than any nonzero change.",
      "refs": ""
    }
  ],
  "script": "import warnings\nimport sys\nimport math\nimport json\nfrom pathlib import Path\n\nimport pandas as pd\nimport numpy as np\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nfrom pydeseq2.dds import DeseqDataSet\nfrom pydeseq2.ds import DeseqStats\nfrom pydeseq2.default_inference import DefaultInference\n\nCOUNTS = \"counts.csv\"\nMETADATA = \"samples.csv\"\n\nEXPECTED = {\n    \"n_samples\": 8,\n    \"n_genes_tested\": 2854,\n    \"n_significant\": 133,\n    \"n_up\": 73,\n    \"n_down\": 60,\n    \"size_factor__control_1\": 0.886627854248,\n    \"size_factor__control_2\": 1.26241510223,\n    \"size_factor__control_3\": 0.932685154622,\n    \"size_factor__control_4\": 1.13362762897,\n    \"size_factor__knockdown_1\": 1.14728149189,\n    \"size_factor__knockdown_2\": 0.732541857301,\n    \"size_factor__knockdown_3\": 1.04549178034,\n    \"size_factor__knockdown_4\": 1.0539813074,\n    \"trend_a0\": 0.02533435,\n    \"trend_a1\": 2.338907,\n    \"prior_disp_var\": 0.25,\n    \"baseMean__GENE1253\": 16186.89507,\n    \"log2FoldChange__GENE1253\": 3.9054175,\n    \"pvalue__GENE1253\": 8.7603854e-144,\n    \"padj__GENE1253\": 2.3057334e-140,\n    \"baseMean__GENE1667\": 502.7724647,\n    \"log2FoldChange__GENE1667\": 2.9812046,\n    \"pvalue__GENE1667\": 1.5706904e-63,\n    \"padj__GENE1667\": 2.0670285e-60,\n    \"baseMean__GENE0757\": 1309.593327,\n    \"log2FoldChange__GENE0757\": 2.4250098,\n    \"pvalue__GENE0757\": 1.6964571e-60,\n    \"padj__GENE0757\": 1.4883584e-57,\n    \"baseMean__GENE0769\": 364.9591817,\n    \"log2FoldChange__GENE0769\": 3.095341,\n    \"pvalue__GENE0769\": 1.3373361e-53,\n    \"padj__GENE0769\": 8.7996716e-51,\n    \"baseMean__GENE0934\": 586.2522152,\n    \"log2FoldChange__GENE0934\": 2.9803354,\n    \"pvalue__GENE0934\": 1.1348048e-50,\n    \"padj__GENE0934\": 5.9736123e-48,\n}\n\n\ndef load_data():\n    counts_path = Path(COUNTS)\n    metadata_path = Path(METADATA)\n    if not counts_path.is_file():\n        sys.exit(f\"Missing counts file: {COUNTS}\")\n    if not metadata_path.is_file():\n        sys.exit(f\"Missing sample sheet file: {METADATA}\")\n    counts_df = pd.read_csv(counts_path, sep=\",\", index_col=0, comment=\"#\")\n    counts_df = counts_df.drop(columns=[], errors=\"ignore\")\n    counts_df = counts_df.T\n    metadata = pd.read_csv(metadata_path, sep=\",\", index_col=0)\n    shared = [s for s in counts_df.index if s in metadata.index]\n    counts_df = counts_df.loc[shared]\n    metadata = metadata.loc[shared]\n    keep = metadata[\"condition\"].notna() & metadata[\"batch\"].notna()\n    metadata = metadata.loc[keep]\n    counts_df = counts_df.loc[metadata.index]\n    metadata[\"batch\"] = metadata[\"batch\"].astype(str)\n    metadata[\"condition\"] = pd.Categorical(\n        metadata[\"condition\"].astype(str), categories=[\"control\", \"knockdown\"]\n    )\n    counts_df = counts_df.loc[:, counts_df.sum(axis=0) >= 10]\n    return counts_df, metadata\n\n\ndef rerun(counts_df, metadata, design, contrast, lfc_null=None):\n    dds2 = DeseqDataSet(\n        counts=counts_df, metadata=metadata, design=design,\n        refit_cooks=True, inference=DefaultInference(n_cpus=1),\n    )\n    dds2.deseq2()\n    extra = {\"lfc_null\": lfc_null, \"alt_hypothesis\": \"greaterAbs\"} if lfc_null is not None else {}\n    ds2 = DeseqStats(\n        dds2, contrast=contrast, alpha=0.05, cooks_filter=True,\n        independent_filter=True, inference=DefaultInference(n_cpus=1), **extra,\n    )\n    ds2.summary()\n    return ds2.results_df\n\n\ndef main():\n    warnings.filterwarnings(\"ignore\")\n    counts_df, metadata = load_data()\n    inference = DefaultInference(n_cpus=1)\n\n    dds = DeseqDataSet(\n        counts=counts_df,\n        metadata=metadata,\n        design=\"~batch + condition\",\n        refit_cooks=True,\n        inference=inference,\n    )\n    dds.deseq2()\n\n    ds = DeseqStats(\n        dds,\n        contrast=[\"condition\", \"knockdown\", \"control\"],\n        alpha=0.05,\n        cooks_filter=True,\n        independent_filter=True,\n        inference=inference,\n    )\n    ds.summary()\n    ds.lfc_shrink(coeff=\"condition[T.knockdown]\")\n\n    res = ds.results_df\n    mismatches = []\n\n    def check(name, got, want, **kw):\n        if not math.isclose(got, want, **kw):\n            mismatches.append(f\"{name}: got {got}, expected {want}\")\n\n    def check_count(name, got, want, tol=2):\n        if abs(got - want) > tol:\n            mismatches.append(f\"{name}: got {got}, expected {want} (tol {tol})\")\n\n    n_samples = dds.n_obs\n    n_genes_tested = len(res)\n    sig_mask = res[\"padj\"] < 0.05\n    n_significant = int(sig_mask.sum())\n    n_up = int((sig_mask & (res[\"log2FoldChange\"] > 0)).sum())\n    n_down = int((sig_mask & (res[\"log2FoldChange\"] < 0)).sum())\n\n    if n_samples != EXPECTED[\"n_samples\"]:\n        mismatches.append(f\"n_samples: got {n_samples}, expected {EXPECTED['n_samples']}\")\n    if n_genes_tested != EXPECTED[\"n_genes_tested\"]:\n        mismatches.append(f\"n_genes_tested: got {n_genes_tested}, expected {EXPECTED['n_genes_tested']}\")\n    check_count(\"n_significant\", n_significant, EXPECTED[\"n_significant\"])\n    check_count(\"n_up\", n_up, EXPECTED[\"n_up\"])\n    check_count(\"n_down\", n_down, EXPECTED[\"n_down\"])\n\n    size_factors = dds.obs[\"size_factors\"]\n    for sample in metadata.index:\n        key = f\"size_factor__{sample}\"\n        if key in EXPECTED:\n            check(key, float(size_factors[sample]), EXPECTED[key], rel_tol=1e-9)\n\n    trend = dds.uns[\"trend_coeffs\"]\n    check(\"trend_a0\", float(trend[\"a0\"]), EXPECTED[\"trend_a0\"], rel_tol=0.1)\n    check(\"trend_a1\", float(trend[\"a1\"]), EXPECTED[\"trend_a1\"], rel_tol=0.1)\n    check(\"prior_disp_var\", float(dds.uns[\"prior_disp_var\"]), EXPECTED[\"prior_disp_var\"], rel_tol=0.1)\n\n    for gene in [\"GENE1253\", \"GENE1667\", \"GENE0757\", \"GENE0769\", \"GENE0934\"]:\n        if gene not in res.index:\n            mismatches.append(f\"{gene}: not found in results\")\n            continue\n        check(f\"baseMean__{gene}\", float(res.loc[gene, \"baseMean\"]), EXPECTED[f\"baseMean__{gene}\"], rel_tol=1e-9)\n        check(f\"log2FoldChange__{gene}\", float(res.loc[gene, \"log2FoldChange\"]), EXPECTED[f\"log2FoldChange__{gene}\"], abs_tol=5e-3)\n        got_p, want_p = float(res.loc[gene, \"pvalue\"]), EXPECTED[f\"pvalue__{gene}\"]\n        check(f\"pvalue__{gene}\", -math.log10(max(got_p, 1e-300)), -math.log10(max(want_p, 1e-300)), rel_tol=0.05, abs_tol=0.05)\n        got_a, want_a = float(res.loc[gene, \"padj\"]), EXPECTED[f\"padj__{gene}\"]\n        check(f\"padj__{gene}\", -math.log10(max(got_a, 1e-300)), -math.log10(max(want_a, 1e-300)), rel_tol=0.05, abs_tol=0.05)\n\n    if mismatches:\n        print(\"EXPECTED mismatches:\")\n        for m in mismatches:\n            print(\" -\", m)\n    else:\n        print(\"No EXPECTED mismatches.\")\n    print(json.dumps({\"n_significant\": n_significant, \"n_up\": n_up, \"n_down\": n_down}))\n\n    # save results + significant genes\n    res.to_csv(\"deseq2_results.csv\")\n    sig = res[res[\"padj\"] < 0.05]\n    sig.to_csv(\"significant_genes.csv\")\n    print(f\"Saved deseq2_results.csv ({len(res)} genes) and significant_genes.csv ({len(sig)} genes).\")\n\n    # top genes by padj\n    top = res.sort_values(\"padj\").head(20)\n    print(\"Top genes by padj:\")\n    print(top[[\"baseMean\", \"log2FoldChange\", \"pvalue\", \"padj\"]])\n\n    # volcano plot\n    plt.figure(figsize=(6, 5))\n    neglog10padj = -np.log10(res[\"padj\"].fillna(1.0).clip(lower=1e-300))\n    plt.scatter(res[\"log2FoldChange\"], neglog10padj, s=5, alpha=0.5)\n    plt.xlabel(\"log2FoldChange\")\n    plt.ylabel(\"-log10(padj)\")\n    plt.title(\"Volcano plot: knockdown vs control\")\n    plt.savefig(\"volcano.png\", dpi=150)\n    plt.close()\n    print(\"Saved volcano.png\")\n\n    # VST PCA for outliers / hidden batch structure\n    dds.vst()\n    vst_counts = np.asarray(dds.layers[\"vst_counts\"])\n    vst_centered = vst_counts - vst_counts.mean(axis=0)\n    u, s, vt = np.linalg.svd(vst_centered, full_matrices=False)\n    pcs = u * s\n    print(\"VST PCA (PC1, PC2) by sample:\")\n    for i, sample in enumerate(metadata.index):\n        cond = metadata[\"condition\"][sample]\n        batch = metadata[\"batch\"][sample]\n        print(f\"  {sample} (condition={cond}, batch={batch}): PC1={pcs[i, 0]:.3f}, PC2={pcs[i, 1]:.3f}\")\n\n    # re-run without the batch term\n    res_nb = rerun(counts_df, metadata, \"~condition\", [\"condition\", \"knockdown\", \"control\"])\n    n_sig_nb = int((res_nb[\"padj\"] < 0.05).sum())\n    print(f\"Without batch term: {n_sig_nb} significant genes (was {n_significant} with ~batch + condition).\")\n\n    # thresholded test lfc_null=1\n    res_th = rerun(counts_df, metadata, \"~batch + condition\", [\"condition\", \"knockdown\", \"control\"], lfc_null=1.0)\n    n_sig_th = int((res_th[\"padj\"] < 0.05).sum())\n    print(f\"Thresholded test (lfc_null=1, greaterAbs): {n_sig_th} significant genes (was {n_significant} at lfc_null=0).\")\n\n\nif __name__ == \"__main__\":\n    main()\n",
  "assumptions": [
    "counts.csv and samples.csv in the working directory are the exact files the page read, unmodified.",
    "pydeseq2 0.5.4 is installed, matching settings.pydeseq2_version.",
    "The values in counts.csv are raw integer read counts, not already normalized or transformed.",
    "The sample and gene identifiers in the two files match exactly between counts.csv and samples.csv."
  ],
  "checks": [
    "The reproduction prints \"No EXPECTED mismatches.\"; a rare mismatch in trend_a0, trend_a1, prior_disp_var or a single gene's pvalue/padj can come from genes whose gene-wise dispersion sits at pydeseq2's lower dispersion bound, where its own fit is decided by rounding.",
    "How n_significant, n_up and n_down (133, 73, 60, M4) move after the no-batch re-run and after the lfc_null=1 re-run.",
    "Whether the VST PCA separates the 8 samples (4 control, 4 knockdown) by condition or by batch (b1, b2) rather than randomly.",
    "That significant_genes.csv has the same row count as the n_significant the script prints."
  ],
  "next_steps": [
    "Run the script and confirm it prints no EXPECTED mismatch against the browser's 36 expected values.",
    "Review the VST PCA printout to see whether the 8 samples separate by condition (4 vs 4) or by batch (b1/b2).",
    "Compare the significant gene count with and without the batch term (currently 133, M4) to judge whether ~batch + condition is warranted.",
    "Send the 133 significant genes (M4) to an enrichment or pathway analysis, since this result is statistical only.",
    "Consider more replicates per condition in a follow-up experiment to raise the residual degrees of freedom above the current 5 (M1)."
  ],
  "prescan_responses": [
    {
      "ref": "F1",
      "verdict": "confirmed",
      "note": "The dispersion prior variance sits exactly at its 0.25 floor (M7); the script's prior_disp_var check reproduces it. It is a low-severity flag and does not change the sound verdict, but genes near the floor are the likely source of any rare rounding-driven mismatch the script reports."
    }
  ]
}