OSCR

Machine learning-based gait classification and genome-wide association identify a QTL for gait type in Colombian paso horses.

Code ↔ Paper

4 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 4 matches
  1. [1] § STAR★Methods › Quantification and statistical analysis › Multi-model GWAS and robustness testing ↔ .opencode/gcta-mcp/server.py, lines 1292–1375 · score 0.80 · Principal Component, population structure, genome wide association, linear model, Mixed Model, Regression
  2. [2] § STAR★Methods › Quantification and statistical analysis › MLMA-LOCO and SAIGE analysis ↔ .opencode/gcta-mcp/server.py, lines 1292–1375 · score 0.73 · Mixed Linear Model, MLMA LOCO, relationship matrix, Chromosome, GCTA, binary
  3. [3] § STAR★Methods › Quantification and statistical analysis › Genotype quality control and filtering ↔ src/Main.cpp, lines 424–466 · score 0.53 · Quality control, allele frequency, QC, MAF
  4. [4] § STAR★Methods › Quantification and statistical analysis › Cross-validation: Comparative frequency analysis ↔ R/SAIGE_fitGLMM_fast.R, lines 240–314 · score 0.51 · Generalized Linear Mixed, fitted, coefficients, Variation, phenotypes, Model

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 2,155 lines · 72 KB · GPL-3.0 · 2 matches

  1. #!/usr/bin/env python3
  2. """
  3. MCP Server for GCTA (Genome-wide Complex Trait Analysis).
  5. Exposes all GCTA functionality through MCP tools so that AI assistants
  6. like Claude, Codex, and opencode can run GCTA analyses via natural language.
  8. Usage:
  9. python server.py
  11. Configuration via environment variables:
  12. GCTA_BINARY_PATH Path to the GCTA binary
  13. (default: ../gcta relative to this script)
  14. GCTA_WORK_DIR Working directory for GCTA data files
  15. (default: parent directory of this script)
  16. GCTA_TIMEOUT Execution timeout in seconds (default: 3600)
  17. """
  19. import os
  20. import sys
  21. import json
  22. import subprocess
  23. import glob
  25. from mcp.server.fastmcp import FastMCP
  28. # ============================================================================
  29. # Configuration
  30. # ============================================================================
  32. _SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
  33. _PKG_ROOT = os.path.dirname(_SCRIPT_DIR)
  35. GCTA_BINARY_PATH = os.environ.get(
  36. "GCTA_BINARY_PATH",
  37. os.path.join(_PKG_ROOT, "gcta"),
  38. )
  39. GCTA_WORK_DIR = os.environ.get(
  40. "GCTA_WORK_DIR",
  41. _PKG_ROOT,
  42. )
  43. GCTA_TIMEOUT = int(os.environ.get("GCTA_TIMEOUT", "3600"))
  46. # ============================================================================
  47. # Execution wrapper
  48. # ============================================================================
  50. class GCTAExecutor:
  51. """Handles execution of GCTA binary via subprocess."""
  53. def __init__(self, binary_path, work_dir, timeout):
  54. self.binary_path = binary_path
  55. self.work_dir = work_dir
  56. self.timeout = timeout
  58. def execute(self, args: str, work_dir: str = None) -> dict:
  59. """Execute GCTA with the given command-line arguments."""
  60. wd = os.path.abspath(work_dir or self.work_dir)
  61. cmd = [self.binary_path] + args.split()
  62. try:
  63. result = subprocess.run(
  64. cmd,
  65. capture_output=True,
  66. text=True,
  67. cwd=wd,
  68. timeout=self.timeout,
  69. )
  70. return {
  71. "exit_code": result.returncode,
  72. "stdout": result.stdout,
  73. "stderr": result.stderr,
  74. "command": " ".join(cmd),
  75. "work_dir": wd,
  76. }
  77. except subprocess.TimeoutExpired:
  78. return {"error": f"GCTA execution timed out after {self.timeout}s"}
  79. except FileNotFoundError:
  80. return {"error": f"GCTA binary not found at {self.binary_path}"}
  81. except Exception as e:
  82. return {"error": str(e)}
  85. # ============================================================================
  86. # Helper functions
  87. # ============================================================================
  89. def _parse_out_prefix(args: str) -> str:
  90. """Parse the --out value from args; return default if absent."""
  91. tokens = args.split()
  92. for i, t in enumerate(tokens):
  93. if t == "--out" and i + 1 < len(tokens):
  94. return tokens[i + 1]
  95. return "gcta"
  98. def _read_log_file(work_dir: str, out_prefix: str) -> str:
  99. """Read the GCTA .log file if it exists."""
  100. log_path = os.path.join(work_dir, out_prefix + ".log")
  101. if os.path.exists(log_path):
  102. try:
  103. with open(log_path, "r", errors="replace") as f:
  104. return f.read()
  105. except Exception:
  106. return ""
  107. return ""
  110. def _list_output_files(work_dir: str, out_prefix: str) -> list:
  111. """List output files matching the out prefix."""
  112. files = []
  113. pattern = os.path.join(work_dir, out_prefix + "*")
  114. for f in sorted(glob.glob(pattern)):
  115. files.append({
  116. "file": os.path.relpath(f, work_dir),
  117. "size_bytes": os.path.getsize(f),
  118. })
  119. return files
  122. def _format_result(exec_result: dict, work_dir: str,
  123. out_prefix: str) -> str:
  124. """Format the execution result as a JSON string."""
  125. if "error" in exec_result:
  126. return json.dumps({"error": exec_result["error"]}, indent=2)
  128. log_content = _read_log_file(work_dir, out_prefix)
  129. output_files = _list_output_files(work_dir, out_prefix)
  131. result = {
  132. "exit_code": exec_result["exit_code"],
  133. "success": exec_result["exit_code"] == 0,
  134. "command": exec_result["command"],
  135. "work_dir": work_dir,
  136. "out_prefix": out_prefix,
  137. "stdout": exec_result["stdout"][-5000:] if exec_result["stdout"] else "",
  138. "stderr": exec_result["stderr"][-5000:] if exec_result["stderr"] else "",
  139. "log_content": log_content[-8000:] if log_content else "",
  140. "output_files": output_files,
  141. }
  142. return json.dumps(result, indent=2, ensure_ascii=False)
  145. def _build_args(base: list, **kwargs) -> str:
  146. """Build a GCTA argument string from a base list and keyword arguments.
  148. ``base`` is a list of required flags + values (e.g.
  149. ``["--bfile", bfile, "--make-grm", "--out", out]``).
  151. Keyword arguments whose value is:
  152. - ``None`` or ``""`` -> skipped
  153. - ``False`` -> skipped (boolean flags are added only when True)
  154. - ``True`` -> added as a bare flag (``--flag``)
  155. - any other value -> added as ``--flag value``
  156. """
  157. parts = list(base)
  158. for key, val in kwargs.items():
  159. if val is None or val == "" or val is False:
  160. continue
  161. flag = "--" + key.replace("_", "-")
  162. if val is True:
  163. parts.append(flag)
  164. else:
  165. parts.extend([flag, str(val)])
  166. return " ".join(parts)
  169. # ============================================================================
  170. # Create MCP server
  171. # ============================================================================
  173. mcp = FastMCP("gcta-mcp")
  175. _executor = GCTAExecutor(
  176. binary_path=GCTA_BINARY_PATH,
  177. work_dir=GCTA_WORK_DIR,
  178. timeout=GCTA_TIMEOUT,
  179. )
  182. # ============================================================================
  183. # Core tools
  184. # ============================================================================
  186. @mcp.tool()
  187. def gcta_run(args: str, work_dir: str = "") -> str:
  188. """Execute GCTA with arbitrary command-line arguments.
  190. This is the most flexible tool — it gives full access to every GCTA
  191. option. Pass everything that would appear after ``gcta64`` on the
  192. command line.
  194. Examples:
  195. args = "--bfile test --make-grm --out test_grm"
  196. args = "--grm test_grm --reml --pheno test.phen --out test_reml"
  197. args = "--bfile test --mlma-loco --pheno test.phen --out test_mlma"
  199. Args:
  200. args: Command-line arguments (everything after the binary name).
  201. work_dir: Working directory for GCTA. Defaults to GCTA_WORK_DIR.
  202. All input/output file paths are relative to this directory.
  204. Returns:
  205. JSON with exit_code, stdout, stderr, log_content, output_files.
  206. """
  207. wd = work_dir or GCTA_WORK_DIR
  208. out_prefix = _parse_out_prefix(args)
  209. result = _executor.execute(args, wd)
  210. return _format_result(result, os.path.abspath(wd), out_prefix)
  213. @mcp.tool()
  214. def gcta_help(topic: str = "all") -> str:
  215. """Get comprehensive documentation about GCTA analyses and options.
  217. Args:
  218. topic: One of:
  219. "all" - Overview of all analyses (default)
  220. "data" - Data input, management, and filtering
  221. "grm" - Genetic Relationship Matrix (GRM) construction
  222. "reml" - REML / GREML variance component estimation
  223. "bivariate_reml" - Bivariate REML analysis
  224. "hereg" - Haseman-Elston regression
  225. "pca" - Principal component analysis
  226. "mlma" - Mixed linear model association (MLMA / MLMA-LOCO)
  227. "cojo" - Conditional and joint analysis (COJO)
  228. "gsmr" - GSMR Mendelian randomization
  229. "mtcojo" - Multi-trait COJO
  230. "fastgwa" - fastGWA genome-wide association
  231. "fastbat" - fastBAT / mBAT gene-based test
  232. "simulation" - Phenotype simulation
  233. "ld" - LD pruning, LD score regression
  234. "fst" - Fst population differentiation
  235. "expression" - Expression data analysis (eRFile, ecojo, make-erm)
  236. "acat" - ACAT gene-based test
  237. "options" - Full list of all command-line options
  239. Returns:
  240. Documentation string for the requested topic.
  241. """
  242. docs = _get_help_docs()
  243. topic = topic.lower().strip()
  244. if topic in docs:
  245. return docs[topic]
  246. elif topic == "all":
  247. return docs["all"]
  248. else:
  249. available = ", ".join(sorted(docs.keys()))
  250. return f"Unknown topic '{topic}'. Available topics: {available}"
  253. @mcp.tool()
  254. def gcta_check() -> str:
  255. """Check if GCTA is properly configured and can be executed.
  257. Verifies:
  258. - GCTA binary exists at the configured path
  259. - Working directory exists and is writable
  260. - Test data files are present
  262. Returns:
  263. JSON with configuration details and check results.
  264. """
  265. info = {
  266. "binary_path": GCTA_BINARY_PATH,
  267. "binary_exists": os.path.exists(GCTA_BINARY_PATH),
  268. "work_dir": GCTA_WORK_DIR,
  269. "work_dir_exists": os.path.exists(GCTA_WORK_DIR),
  270. "work_dir_writable": os.access(GCTA_WORK_DIR, os.W_OK) if os.path.exists(GCTA_WORK_DIR) else False,
  271. "timeout_seconds": GCTA_TIMEOUT,
  272. }
  274. checks = []
  276. if not info["binary_exists"]:
  277. checks.append("FAIL: GCTA binary not found. Set GCTA_BINARY_PATH env var.")
  278. else:
  279. checks.append("OK: GCTA binary found.")
  281. if not info["work_dir_exists"]:
  282. checks.append("FAIL: Working directory not found. Set GCTA_WORK_DIR env var.")
  283. else:
  284. checks.append("OK: Working directory exists.")
  285. if info["work_dir_writable"]:
  286. checks.append("OK: Working directory is writable.")
  287. else:
  288. checks.append("WARN: Working directory is not writable.")
  290. test_data_dir = GCTA_WORK_DIR
  291. test_files = ["test.bed", "test.bim", "test.fam", "test.phen"]
  292. missing = [f for f in test_files if not os.path.exists(os.path.join(test_data_dir, f))]
  293. if not missing:
  294. checks.append("OK: Test data files found in working directory.")
  295. else:
  296. checks.append(f"WARN: Missing test data files: {', '.join(missing)}")
  298. info["checks"] = checks
  299. return json.dumps(info, indent=2, ensure_ascii=False)
  302. # ============================================================================
  303. # GRM construction tools
  304. # ============================================================================
  306. @mcp.tool()
  307. def gcta_make_grm(
  308. bfile: str,
  309. out: str = "gcta_grm",
  310. maf: float = 0.0,
  311. thread_num: int = 1,
  312. autosome_num: int = 22,
  313. chr_filter: int = 0,
  314. autosome: bool = False,
  315. make_grm_alg: int = 0,
  316. grm_adj: float = -1.0,
  317. grm_cutoff: float = -1.0,
  318. dosage_compen: int = -1,
  319. dominance: bool = False,
  320. make_grm_xchr: bool = False,
  321. make_grm_inbred: bool = False,
  322. save_ram: bool = False,
  323. extract: str = "",
  324. exclude: str = "",
  325. work_dir: str = "",
  326. ) -> str:
  327. """Construct a Genetic Relationship Matrix (GRM) from PLINK bed data.
  329. The GRM estimates pairwise genomic relatedness between individuals
  330. based on genome-wide SNPs. This is the foundational step for most
  331. downstream GCTA analyses (REML, PCA, MLMA, etc.).
  333. Args:
  334. bfile: PLINK binary file prefix (e.g. "test" for test.bed/.bim/.fam).
  335. out: Output file prefix. Produces {out}.grm.bin, {out}.grm.N.bin,
  336. {out}.grm.id.
  337. maf: Minor allele frequency filter (0 = no filter, typical: 0.01).
  338. thread_num: Number of CPU threads.
  339. autosome_num: Number of autosomes (default 22 for human).
  340. chr_filter: Restrict to a specific chromosome (0 = no restriction).
  341. autosome: If True, restrict to autosomes 1..autosome_num.
  342. make_grm_alg: Algorithm: 0 = VanRaden (default), 1 = alternative.
  343. grm_adj: GRM adjustment factor (-1 = not used).
  344. grm_cutoff: Remove pairs with GRM > cutoff (-1 = no cutoff).
  345. dosage_compen: Dosage compensation: 0 = none, 1 = X chromosome.
  346. dominance: If True, calculate dominance GRM.
  347. make_grm_xchr: If True, calculate X-chromosome GRM.
  348. make_grm_inbred: If True, calculate inbreeding GRM.
  349. save_ram: If True, use memory-saving mode.
  350. extract: File of SNP IDs to include.
  351. exclude: File of SNP IDs to exclude.
  352. work_dir: Working directory for GCTA.
  354. Returns:
  355. JSON with execution results, log content, and output file list.
  356. """
  357. base = ["--bfile", bfile, "--autosome-num", str(autosome_num),
  358. "--thread-num", str(thread_num)]
  360. if dominance:
  361. base.append("--make-grm-d")
  362. elif make_grm_xchr:
  363. base.append("--make-grm-xchr")
  364. elif make_grm_inbred:
  365. base.append("--make-grm-inbred")
  366. else:
  367. base.append("--make-grm")
  369. base.extend(["--out", out])
  371. args = _build_args(base, maf=maf, chr=chr_filter if chr_filter > 0 else None,
  372. autosome=autosome, make_grm_alg=make_grm_alg,
  373. grm_adj=grm_adj if grm_adj >= 0 else None,
  374. grm_cutoff=grm_cutoff if grm_cutoff > 0 else None,
  375. dc=dosage_compen if dosage_compen >= 0 else None,
  376. save_ram=save_ram, extract=extract or None,
  377. exclude=exclude or None)
  378. return gcta_run(args, work_dir)
  381. # ============================================================================
  382. # REML / GREML tools
  383. # ============================================================================
  385. @mcp.tool()
  386. def gcta_reml(
  387. grm: str,
  388. pheno: str,
  389. out: str = "gcta_reml",
  390. qcovar: str = "",
  391. covar: str = "",
  392. mpheno: int = 0,
  393. thread_num: int = 1,
  394. reml_priors: str = "",
  395. reml_priors_var: str = "",
  396. reml_maxit: int = 100,
  397. prevalence: float = -1.0,
  398. reml_no_constrain: bool = False,
  399. reml_no_lrt: bool = False,
  400. reml_pred_rand: bool = False,
  401. reml_est_fix: bool = False,
  402. grm_cutoff: float = -1.0,
  403. gxe: str = "",
  404. work_dir: str = "",
  405. ) -> str:
  406. """Run REML analysis to estimate variance components (GREML).
  408. Estimates the proportion of phenotypic variance explained by all
  409. genome-wide SNPs (SNP-based heritability, h2_SNP) using Restricted
  410. Maximum Likelihood (REML).
  412. Args:
  413. grm: GRM file prefix (from gcta_make_grm).
  414. pheno: Phenotype file path.
  415. out: Output file prefix. Produces {out}.hsq.
  416. qcovar: Quantitative covariate file.
  417. covar: Categorical covariate file.
  418. mpheno: Phenotype column number (0 = first column).
  419. thread_num: Number of CPU threads.
  420. reml_priors: Prior values of variance explained (space-separated).
  421. reml_priors_var: Prior values of variance components (space-separated).
  422. reml_maxit: Maximum REML iterations (default 100).
  423. prevalence: Disease prevalence for liability conversion (-1 = not used).
  424. reml_no_constrain: If True, do not constrain variance components.
  425. reml_no_lrt: If True, skip likelihood ratio test.
  426. reml_pred_rand: If True, predict random effects (BLUP).
  427. reml_est_fix: If True, estimate fixed effects.
  428. grm_cutoff: GRM cutoff for related pairs (-1 = no cutoff).
  429. gxe: GxE interaction file.
  430. work_dir: Working directory for GCTA.
  432. Returns:
  433. JSON with execution results, log content, and output file list.
  434. """
  435. base = ["--grm", grm, "--pheno", pheno,
  436. "--thread-num", str(thread_num), "--reml",
  437. "--reml-maxit", str(reml_maxit), "--out", out]
  439. args = _build_args(base, qcovar=qcovar or None, covar=covar or None,
  440. mpheno=mpheno or None,
  441. reml_priors=reml_priors or None,
  442. reml_priors_var=reml_priors_var or None,
  443. prevalence=prevalence if prevalence > 0 else None,
  444. reml_no_constrain=reml_no_constrain,
  445. reml_no_lrt=reml_no_lrt,
  446. reml_pred_rand=reml_pred_rand,
  447. reml_est_fix=reml_est_fix,
  448. grm_cutoff=grm_cutoff if grm_cutoff > 0 else None,
  449. gxe=gxe or None)
  450. return gcta_run(args, work_dir)
  453. @mcp.tool()
  454. def gcta_bivariate_reml(
  455. grm: str,
  456. pheno: str,
  457. out: str = "gcta_bireml",
  458. qcovar: str = "",
  459. covar: str = "",
  460. mphen: int = 1,
  461. mphen2: int = 2,
  462. thread_num: int = 1,
  463. reml_maxit: int = 100,
  464. reml_bivar_prevalence: str = "",
  465. reml_bivar_nocove: bool = False,
  466. reml_bivar_no_constrain: bool = False,
  467. grm_cutoff: float = -1.0,
  468. work_dir: str = "",
  469. ) -> str:
  470. """Run bivariate REML analysis to estimate genetic correlation (rg).
  472. Estimates the genetic correlation between two traits using bivariate
  473. GREML. Requires a phenotype file with at least two columns.
  475. Args:
  476. grm: GRM file prefix.
  477. pheno: Phenotype file (with at least 2 trait columns).
  478. out: Output prefix. Produces {out}.hsq.
  479. qcovar: Quantitative covariate file.
  480. covar: Categorical covariate file.
  481. mphen: First phenotype column number (default 1).
  482. mphen2: Second phenotype column number (default 2).
  483. thread_num: Number of CPU threads.
  484. reml_maxit: Maximum REML iterations.
  485. reml_bivar_prevalence: Two prevalence values (space-separated).
  486. reml_bivar_nocove: If True, ignore residual covariance.
  487. reml_bivar_no_constrain: If True, no constrain on variance components.
  488. grm_cutoff: GRM cutoff (-1 = no cutoff).
  489. work_dir: Working directory for GCTA.
  491. Returns:
  492. JSON with execution results, log content, and output file list.
  493. """
  494. base = ["--grm", grm, "--pheno", pheno,
  495. "--thread-num", str(thread_num),
  496. "--reml-bivar", str(mphen), str(mphen2),
  497. "--reml-maxit", str(reml_maxit), "--out", out]
  499. args = _build_args(base, qcovar=qcovar or None, covar=covar or None,
  500. reml_bivar_prevalence=reml_bivar_prevalence or None,
  501. reml_bivar_nocove=reml_bivar_nocove,
  502. reml_bivar_no_constrain=reml_bivar_no_constrain,
  503. grm_cutoff=grm_cutoff if grm_cutoff > 0 else None)
  504. return gcta_run(args, work_dir)
  507. @mcp.tool()
  508. def gcta_hereg(
  509. grm: str,
  510. pheno: str,
  511. out: str = "gcta_hereg",
  512. qcovar: str = "",
  513. covar: str = "",
  514. mpheno: int = 0,
  515. thread_num: int = 1,
  516. grm_cutoff: float = -1.0,
  517. work_dir: str = "",
  518. ) -> str:
  519. """Run Haseman-Elston regression to estimate variance components.
  521. HE regression is a method-of-moments estimator that is faster but
  522. less precise than REML. Useful as a sanity check or when REML
  523. fails to converge.
  525. Args:
  526. grm: GRM file prefix.
  527. pheno: Phenotype file path.
  528. out: Output prefix.
  529. qcovar: Quantitative covariate file.
  530. covar: Categorical covariate file.
  531. mpheno: Phenotype column number (0 = first).
  532. thread_num: Number of CPU threads.
  533. grm_cutoff: GRM cutoff (-1 = no cutoff).
  534. work_dir: Working directory for GCTA.
  536. Returns:
  537. JSON with execution results, log content, and output file list.
  538. """
  539. base = ["--grm", grm, "--pheno", pheno,
  540. "--thread-num", str(thread_num), "--HEreg",
  541. "--out", out]
  543. args = _build_args(base, qcovar=qcovar or None, covar=covar or None,
  544. mpheno=mpheno or None,
  545. grm_cutoff=grm_cutoff if grm_cutoff > 0 else None)
  546. return gcta_run(args, work_dir)
  549. # ============================================================================
  550. # PCA tools
  551. # ============================================================================
  553. @mcp.tool()
  554. def gcta_pca(
  555. grm: str,
  556. out: str = "gcta_pca",
  557. pc_num: int = 20,
  558. thread_num: int = 1,
  559. work_dir: str = "",
  560. ) -> str:
  561. """Perform Principal Component Analysis (PCA) from a GRM.
  563. Computes principal components (PCs) from the genetic relationship
  564. matrix to capture population structure. Output eigenvalues and
  565. eigenvectors are saved in {out}.eigenval and {out}.eigenvec.
  567. Args:
  568. grm: GRM file prefix (from gcta_make_grm).
  569. out: Output file prefix.
  570. pc_num: Number of principal components to output (default 20).
  571. thread_num: Number of CPU threads.
  572. work_dir: Working directory for GCTA.
  574. Returns:
  575. JSON with execution results, log content, and output file list.
  576. """
  577. args = (f"--grm {grm} --pca {pc_num} "
  578. f"--thread-num {thread_num} --out {out}")
  579. return gcta_run(args, work_dir)
  582. # ============================================================================
  583. # Association analysis tools
  584. # ============================================================================
  586. @mcp.tool()
  587. def gcta_mlma(
  588. bfile: str,
  589. grm: str,
  590. pheno: str,
  591. out: str = "gcta_mlma",
  592. qcovar: str = "",
  593. covar: str = "",
  594. mpheno: int = 0,
  595. thread_num: int = 1,
  596. reml_maxit: int = 100,
  597. loco: bool = False,
  598. maf: float = 0.0,
  599. autosome_num: int = 22,
  600. work_dir: str = "",
  601. ) -> str:
  602. """Run Mixed Linear Model Association (MLMA) analysis.
  604. Performs genome-wide association analysis using a linear mixed model
  605. that accounts for population structure via the GRM.
  607. Set ``loco=True`` to use the Leave-One-Chromosome-Out (LOCO) method,
  608. which is more computationally efficient for large datasets.
  610. Args:
  611. bfile: PLINK binary file prefix.
  612. grm: GRM file prefix.
  613. pheno: Phenotype file path.
  614. out: Output file prefix. Produces {out}.mlma (or .loco.mlma).
  615. qcovar: Quantitative covariate file.
  616. covar: Categorical covariate file.
  617. mpheno: Phenotype column number (0 = first).
  618. thread_num: Number of CPU threads.
  619. reml_maxit: Maximum REML iterations.
  620. loco: If True, use MLMA-LOCO method.
  621. maf: Minor allele frequency filter (0 = no filter).
  622. autosome_num: Number of autosomes.
  623. work_dir: Working directory for GCTA.
  625. Returns:
  626. JSON with execution results, log content, and output file list.
  627. """
  628. if loco:
  629. mlma_flag = "--mlma-loco"
  630. else:
  631. mlma_flag = "--mlma"
  633. base = ["--bfile", bfile, "--grm", grm, "--pheno", pheno,
  634. "--autosome-num", str(autosome_num),
  635. "--thread-num", str(thread_num),
  636. mlma_flag, "--reml-maxit", str(reml_maxit), "--out", out]
  638. args = _build_args(base, qcovar=qcovar or None, covar=covar or None,
  639. mpheno=mpheno or None, maf=maf if maf > 0 else None)
  640. return gcta_run(args, work_dir)
  643. @mcp.tool()
  644. def gcta_cojo(
  645. bfile: str,
  646. cojo_file: str,
  647. out: str = "gcta_cojo",
  648. method: str = "slct",
  649. cojo_p: float = 5e-8,
  650. cojo_wind: int = 10000,
  651. cojo_collinear: float = 0.9,
  652. cojo_cond_snplist: str = "",
  653. cojo_gc: float = -1.0,
  654. cojo_sblup_fac: float = -1.0,
  655. cojo_top_snps: int = -1,
  656. thread_num: int = 1,
  657. maf: float = 0.0,
  658. work_dir: str = "",
  659. ) -> str:
  660. """Run Conditional and Joint (COJO) analysis of GWAS summary data.
  662. COJO performs stepwise selection, conditional analysis, or joint
  663. analysis of GWAS summary statistics, using an LD reference panel
  664. from PLINK bed data.
  666. Methods:
  667. "slct" - Stepwise selection (backward + forward)
  668. "forward" - Forward-only selection
  669. "backward" - Backward-only selection
  670. "cond" - Conditional analysis (condition on given SNPs)
  671. "joint" - Joint analysis (all SNPs in one model)
  672. "sblup" - SBLUP (summary-data BLUP)
  674. Args:
  675. bfile: PLINK binary file prefix (LD reference panel).
  676. cojo_file: GWAS summary data file (.ma format).
  677. out: Output file prefix.
  678. method: COJO method (see above).
  679. cojo_p: P-value threshold for selection (default 5e-8).
  680. cojo_wind: Window size in Kb for LD calculation (default 10000).
  681. cojo_collinear: Collinearity threshold (default 0.9).
  682. cojo_cond_snplist: SNP list file for conditional analysis.
  683. cojo_gc: Genomic control inflation factor (-1 = not used).
  684. cojo_sblup_fac: SBLUP factor (for sblup method).
  685. cojo_top_snps: Number of top SNPs to select (-1 = no limit).
  686. thread_num: Number of CPU threads.
  687. maf: Minor allele frequency filter.
  688. work_dir: Working directory for GCTA.
  690. Returns:
  691. JSON with execution results, log content, and output file list.
  692. """
  693. method = method.lower().strip()
  694. method_map = {
  695. "slct": "--cojo-slct",
  696. "stepwise": "--cojo-slct",
  697. "forward": "--cojo-forward",
  698. "backward": "--cojo-backward",
  699. "joint": "--cojo-joint",
  700. "sblup": "--cojo-sblup",
  701. }
  702. if method not in method_map and method != "cond":
  703. return json.dumps({
  704. "error": f"Unknown method '{method}'. Use: slct, forward, backward, cond, joint, sblup"
  705. }, indent=2)
  707. base = ["--bfile", bfile, "--cojo-file", cojo_file,
  708. "--thread-num", str(thread_num), "--out", out]
  710. if method == "cond":
  711. base.extend(["--cojo-cond", cojo_cond_snplist])
  712. elif method == "sblup":
  713. base.extend([method_map[method], str(cojo_sblup_fac)])
  714. else:
  715. base.append(method_map[method])
  717. args = _build_args(base, cojo_p=cojo_p, cojo_wind=cojo_wind,
  718. cojo_collinear=cojo_collinear,
  719. cojo_gc=cojo_gc if cojo_gc > 0 else None,
  720. cojo_top_snps=cojo_top_snps if cojo_top_snps > 0 else None,
  721. maf=maf if maf > 0 else None)
  722. return gcta_run(args, work_dir)
  725. @mcp.tool()
  726. def gcta_gsmr(
  727. bfile: str,
  728. gsmr_file: str,
  729. out: str = "gcta_gsmr",
  730. gsmr_direction: int = 0,
  731. gsmr2_beta: bool = False,
  732. thread_num: int = 1,
  733. maf: float = 0.0,
  734. work_dir: str = "",
  735. ) -> str:
  736. """Run GSMR (Generalized Summary-data Mendelian Randomization) analysis.
  738. GSMR estimates the causal effect of an exposure on an outcome using
  739. GWAS summary data and an LD reference panel.
  741. Args:
  742. bfile: PLINK binary file prefix (LD reference).
  743. gsmr_file: File listing exposure and outcome GWAS summary files.
  744. out: Output file prefix.
  745. gsmr_direction: 0 = forward-GSMR, 1 = reverse-GSMR, 2 = bi-GSMR.
  746. gsmr2_beta: If True, use GSMR2-beta version (multi-SNP HEIDI).
  747. thread_num: Number of CPU threads.
  748. maf: Minor allele frequency filter.
  749. work_dir: Working directory for GCTA.
  751. Returns:
  752. JSON with execution results, log content, and output file list.
  753. """
  754. base = ["--bfile", bfile, "--gsmr-file", gsmr_file,
  755. "--gsmr-direction", str(gsmr_direction),
  756. "--thread-num", str(thread_num), "--out", out]
  758. args = _build_args(base, gsmr2_beta=gsmr2_beta,
  759. maf=maf if maf > 0 else None)
  760. return gcta_run(args, work_dir)
  763. @mcp.tool()
  764. def gcta_mtcojo(
  765. bfile: str,
  766. mtcojo_file: str,
  767. mtcojo_bxy: str,
  768. out: str = "gcta_mtcojo",
  769. thread_num: int = 1,
  770. maf: float = 0.0,
  771. work_dir: str = "",
  772. ) -> str:
  773. """Run mtCOJO (Multi-trait COJO) analysis.
  775. mtCOJO performs conditional analysis across multiple traits using
  776. GWAS summary data, adjusting for the genetic correlation between
  777. traits.
  779. Args:
  780. bfile: PLINK binary file prefix (LD reference).
  781. mtcojo_file: File listing trait GWAS summary files.
  782. mtcojo_bxy: File with causal effect estimates between traits.
  783. out: Output file prefix.
  784. thread_num: Number of CPU threads.
  785. maf: Minor allele frequency filter.
  786. work_dir: Working directory for GCTA.
  788. Returns:
  789. JSON with execution results, log content, and output file list.
  790. """
  791. base = ["--bfile", bfile, "--mtcojo-file", mtcojo_file,
  792. "--mtcojo-bxy", mtcojo_bxy,
  793. "--thread-num", str(thread_num), "--out", out]
  795. args = _build_args(base, maf=maf if maf > 0 else None)
  796. return gcta_run(args, work_dir)
  799. @mcp.tool()
  800. def gcta_fastgwa(
  801. bfile: str,
  802. pheno: str,
  803. out: str = "gcta_fastgwa",
  804. qcovar: str = "",
  805. covar: str = "",
  806. grm: str = "",
  807. method: str = "fastGWA-mlm",
  808. thread_num: int = 1,
  809. maf: float = 0.0,
  810. autosome_num: int = 22,
  811. work_dir: str = "",
  812. ) -> str:
  813. """Run fastGWA genome-wide association analysis.
  815. fastGWA is an efficient mixed-model association method that uses
  816. sparse GRM techniques for improved speed on large datasets.
  818. Methods:
  819. "fastGWA" - fastGWA (requires --grm)
  820. "fastGWA-mlm" - fastGWA-MLM (linear mixed model)
  821. "fastGWA-mlm-exact" - fastGWA-MLM exact mode
  822. "fastGWA-lr" - fastGWA linear regression (no GRM)
  824. Args:
  825. bfile: PLINK binary file prefix.
  826. pheno: Phenotype file path.
  827. out: Output file prefix.
  828. qcovar: Quantitative covariate file.
  829. covar: Categorical covariate file.
  830. grm: GRM file prefix (for fastGWA with GRM).
  831. method: fastGWA method (see above).
  832. thread_num: Number of CPU threads.
  833. maf: Minor allele frequency filter.
  834. autosome_num: Number of autosomes.
  835. work_dir: Working directory for GCTA.
  837. Returns:
  838. JSON with execution results, log content, and output file list.
  839. """
  840. method = method.lower().strip()
  841. method_map = {
  842. "fastgwa": "--fastGWA",
  843. "fastgwa-mlm": "--fastGWA-mlm",
  844. "fastgwa-mlm-exact": "--fastGWA-mlm-exact",
  845. "fastgwa-lr": "--fastGWA-lr",
  846. }
  847. if method not in method_map:
  848. return json.dumps({
  849. "error": f"Unknown method '{method}'. Use: fastGWA, fastGWA-mlm, fastGWA-mlm-exact, fastGWA-lr"
  850. }, indent=2)
  852. base = ["--bfile", bfile, "--pheno", pheno,
  853. "--autosome-num", str(autosome_num),
  854. "--thread-num", str(thread_num),
  855. method_map[method], "--out", out]
  857. args = _build_args(base, qcovar=qcovar or None, covar=covar or None,
  858. grm=grm or None, maf=maf if maf > 0 else None)
  859. return gcta_run(args, work_dir)
  862. @mcp.tool()
  863. def gcta_fastbat(
  864. bfile: str,
  865. fastbat_file: str,
  866. out: str = "gcta_fastbat",
  867. fastbat_gene_list: str = "",
  868. fastbat_set_list: str = "",
  869. fastbat_wind: int = 50,
  870. fastbat_ld_cutoff: float = 0.9,
  871. thread_num: int = 1,
  872. maf: float = 0.0,
  873. work_dir: str = "",
  874. ) -> str:
  875. """Run fastBAT gene-based association test.
  877. fastBAT performs gene-based association tests by combining
  878. SNP-level p-values within genes, accounting for LD structure.
  880. Args:
  881. bfile: PLINK binary file prefix (LD reference).
  882. fastbat_file: GWAS summary data file (.ma format).
  883. out: Output file prefix.
  884. fastbat_gene_list: Gene annotation file (gene -> SNPs mapping).
  885. fastbat_set_list: SNP set file (predefined sets).
  886. fastbat_wind: Window size in Kb for gene boundaries (default 50).
  887. fastbat_ld_cutoff: LD r2 cutoff for removing correlated SNPs (default 0.9).
  888. thread_num: Number of CPU threads.
  889. maf: Minor allele frequency filter.
  890. work_dir: Working directory for GCTA.
  892. Returns:
  893. JSON with execution results, log content, and output file list.
  894. """
  895. base = ["--bfile", bfile, "--fastBAT", fastbat_file,
  896. "--thread-num", str(thread_num), "--out", out]
  898. args = _build_args(base, fastBAT_gene_list=fastbat_gene_list or None,
  899. fastBAT_set_list=fastbat_set_list or None,
  900. fastBAT_wind=fastbat_wind,
  901. fastBAT_ld_cutoff=fastbat_ld_cutoff,
  902. maf=maf if maf > 0 else None)
  903. return gcta_run(args, work_dir)
  906. # ============================================================================
  907. # Simulation tools
  908. # ============================================================================
  910. @mcp.tool()
  911. def gcta_simu_qt(
  912. bfile: str,
  913. out: str = "gcta_simu_qt",
  914. simu_hsq: float = 0.1,
  915. simu_rep: int = 1,
  916. simu_causal_loci: str = "",
  917. simu_seed: float = -1.0,
  918. thread_num: int = 1,
  919. maf: float = 0.0,
  920. work_dir: str = "",
  921. ) -> str:
  922. """Simulate quantitative trait phenotypes based on real genotype data.
  924. Uses the genotype data from PLINK bed files to simulate quantitative
  925. traits with a specified heritability. Useful for power analysis and
  926. method evaluation.
  928. Args:
  929. bfile: PLINK binary file prefix.
  930. out: Output file prefix.
  931. simu_hsq: Simulated heritability (default 0.1).
  932. simu_rep: Number of simulation repetitions (default 1).
  933. simu_causal_loci: File listing causal loci.
  934. simu_seed: Random seed for simulation (-1 = auto).
  935. thread_num: Number of CPU threads.
  936. maf: Minor allele frequency filter.
  937. work_dir: Working directory for GCTA.
  939. Returns:
  940. JSON with execution results, log content, and output file list.
  941. """
  942. base = ["--bfile", bfile, "--simu-qt",
  943. "--simu-hsq", str(simu_hsq),
  944. "--simu-rep", str(simu_rep),
  945. "--thread-num", str(thread_num), "--out", out]
  947. args = _build_args(base, simu_causal_loci=simu_causal_loci or None,
  948. simu_seed=simu_seed if simu_seed > 0 else None,
  949. maf=maf if maf > 0 else None)
  950. return gcta_run(args, work_dir)
  953. @mcp.tool()
  954. def gcta_simu_cc(
  955. bfile: str,
  956. simu_case_num: int,
  957. simu_control_num: int,
  958. out: str = "gcta_simu_cc",
  959. simu_hsq: float = 0.1,
  960. simu_k: float = 0.1,
  961. simu_rep: int = 1,
  962. simu_causal_loci: str = "",
  963. simu_seed: float = -1.0,
  964. thread_num: int = 1,
  965. maf: float = 0.0,
  966. work_dir: str = "",
  967. ) -> str:
  968. """Simulate case-control phenotypes based on real genotype data.
  970. Uses genotype data from PLINK bed files to simulate case-control
  971. phenotypes with specified heritability and disease prevalence.
  973. Args:
  974. bfile: PLINK binary file prefix.
  975. simu_case_num: Number of cases to simulate.
  976. simu_control_num: Number of controls to simulate.
  977. out: Output file prefix.
  978. simu_hsq: Simulated heritability (default 0.1).
  979. simu_k: Simulated disease prevalence (default 0.1).
  980. simu_rep: Number of simulation repetitions (default 1).
  981. simu_causal_loci: File listing causal loci.
  982. simu_seed: Random seed (-1 = auto).
  983. thread_num: Number of CPU threads.
  984. maf: Minor allele frequency filter.
  985. work_dir: Working directory for GCTA.
  987. Returns:
  988. JSON with execution results, log content, and output file list.
  989. """
  990. base = ["--bfile", bfile, "--simu-cc",
  991. str(simu_case_num), str(simu_control_num),
  992. "--simu-hsq", str(simu_hsq),
  993. "--simu-k", str(simu_k),
  994. "--simu-rep", str(simu_rep),
  995. "--thread-num", str(thread_num), "--out", out]
  997. args = _build_args(base, simu_causal_loci=simu_causal_loci or None,
  998. simu_seed=simu_seed if simu_seed > 0 else None,
  999. maf=maf if maf > 0 else None)
  1000. return gcta_run(args, work_dir)
  1003. # ============================================================================
  1004. # Population genetics tools
  1005. # ============================================================================
  1007. @mcp.tool()
  1008. def gcta_fst(
  1009. bfile: str,
  1010. sub_popu: str,
  1011. out: str = "gcta_fst",
  1012. thread_num: int = 1,
  1013. maf: float = 0.0,
  1014. work_dir: str = "",
  1015. ) -> str:
  1016. """Calculate Fst (fixation index) for population differentiation analysis.
  1018. Computes pairwise Fst between subpopulations defined in the
  1019. sub-population file.
  1021. Args:
  1022. bfile: PLINK binary file prefix.
  1023. sub_popu: Subpopulation assignment file (FID, IID, pop ID).
  1024. out: Output file prefix.
  1025. thread_num: Number of CPU threads.
  1026. maf: Minor allele frequency filter.
  1027. work_dir: Working directory for GCTA.
  1029. Returns:
  1030. JSON with execution results, log content, and output file list.
  1031. """
  1032. args = (f"--bfile {bfile} --fst --sub-popu {sub_popu} "
  1033. f"--thread-num {thread_num} --out {out}")
  1034. if maf > 0:
  1035. args += f" --maf {maf}"
  1036. return gcta_run(args, work_dir)
  1039. # ============================================================================
  1040. # Data management tools
  1041. # ============================================================================
  1043. @mcp.tool()
  1044. def gcta_make_bed(
  1045. bfile: str,
  1046. out: str = "gcta_bed",
  1047. maf: float = 0.0,
  1048. thread_num: int = 1,
  1049. extract: str = "",
  1050. exclude: str = "",
  1051. chr_filter: int = 0,
  1052. autosome: bool = False,
  1053. autosome_num: int = 22,
  1054. work_dir: str = "",
  1055. ) -> str:
  1056. """Convert / filter genotype data and output PLINK bed format.
  1058. Performs data management operations: SNP/individual filtering,
  1059. chromosome extraction, and output in PLINK binary format.
  1061. Args:
  1062. bfile: PLINK binary file prefix (input).
  1063. out: Output file prefix.
  1064. maf: Minor allele frequency filter (0 = no filter).
  1065. thread_num: Number of CPU threads.
  1066. extract: File of SNP IDs to include.
  1067. exclude: File of SNP IDs to exclude.
  1068. chr_filter: Restrict to a specific chromosome (0 = no restriction).
  1069. autosome: If True, restrict to autosomes.
  1070. autosome_num: Number of autosomes.
  1071. work_dir: Working directory for GCTA.
  1073. Returns:
  1074. JSON with execution results, log content, and output file list.
  1075. """
  1076. base = ["--bfile", bfile, "--make-bed",
  1077. "--autosome-num", str(autosome_num),
  1078. "--thread-num", str(thread_num), "--out", out]
  1080. args = _build_args(base, maf=maf if maf > 0 else None,
  1081. extract=extract or None, exclude=exclude or None,
  1082. chr=chr_filter if chr_filter > 0 else None,
  1083. autosome=autosome)
  1084. return gcta_run(args, work_dir)
  1087. @mcp.tool()
  1088. def gcta_ld_pruning(
  1089. bfile: str,
  1090. ld_pruning_rsq: float = 0.1,
  1091. out: str = "gcta_ld_prune",
  1092. ld_wind: float = 10000.0,
  1093. thread_num: int = 1,
  1094. maf: float = 0.0,
  1095. work_dir: str = "",
  1096. ) -> str:
  1097. """Perform LD pruning on genotype data.
  1099. Removes SNPs in high LD with each other, retaining a set of
  1100. approximately independent SNPs. Useful for reducing redundancy
  1101. before downstream analyses.
  1103. Args:
  1104. bfile: PLINK binary file prefix.
  1105. ld_pruning_rsq: LD r2 threshold for pruning (default 0.1).
  1106. out: Output file prefix.
  1107. ld_wind: LD window size in Kb (default 10000).
  1108. thread_num: Number of CPU threads.
  1109. maf: Minor allele frequency filter.
  1110. work_dir: Working directory for GCTA.
  1112. Returns:
  1113. JSON with execution results, log content, and output file list.
  1114. """
  1115. args = (f"--bfile {bfile} --ld-pruning {ld_pruning_rsq} "
  1116. f"--ld-wind {ld_wind} "
  1117. f"--thread-num {thread_num} --out {out}")
  1118. if maf > 0:
  1119. args += f" --maf {maf}"
  1120. return gcta_run(args, work_dir)
  1123. @mcp.tool()
  1124. def gcta_ld_score(
  1125. bfile: str,
  1126. out: str = "gcta_ld_score",
  1127. ld_wind: float = 10000.0,
  1128. thread_num: int = 1,
  1129. maf: float = 0.0,
  1130. work_dir: str = "",
  1131. ) -> str:
  1132. """Calculate LD scores for SNPs in genotype data.
  1134. Computes LD scores (sum of LD r2 with neighboring SNPs) for each
  1135. SNP. Used as input for LD score regression analyses.
  1137. Args:
  1138. bfile: PLINK binary file prefix.
  1139. out: Output file prefix.
  1140. ld_wind: LD window size in Kb (default 10000).
  1141. thread_num: Number of CPU threads.
  1142. maf: Minor allele frequency filter.
  1143. work_dir: Working directory for GCTA.
  1145. Returns:
  1146. JSON with execution results, log content, and output file list.
  1147. """
  1148. args = (f"--bfile {bfile} --ld-score --ld-wind {ld_wind} "
  1149. f"--thread-num {thread_num} --out {out}")
  1150. if maf > 0:
  1151. args += f" --maf {maf}"
  1152. return gcta_run(args, work_dir)
  1155. # ============================================================================
  1156. # ACAT tools
  1157. # ============================================================================
  1159. @mcp.tool()
  1160. def gcta_acat(
  1161. gene_list: str,
  1162. snp_list: str,
  1163. out: str = "acat_res.csv",
  1164. max_maf: float = 0.01,
  1165. min_mac: int = 20,
  1166. wind: int = 0,
  1167. work_dir: str = "",
  1168. ) -> str:
  1169. """Run ACAT (Aggregated Cauchy Association Test) gene-based test.
  1171. ACAT combines p-values from multiple SNPs within a gene using a
  1172. Cauchy combination method, providing a gene-level association test.
  1174. Args:
  1175. gene_list: Gene list file (gene -> chromosome, start, end).
  1176. snp_list: SNP list file with p-values.
  1177. out: Output file (default acat_res.csv).
  1178. max_maf: Maximum MAF for included SNPs (default 0.01).
  1179. min_mac: Minimum minor allele count (default 20).
  1180. wind: Extension length for gene boundaries (default 0).
  1181. work_dir: Working directory for GCTA.
  1183. Returns:
  1184. JSON with execution results, log content, and output file list.
  1185. """
  1186. args = (f"--acat --gene-list {gene_list} --snp-list {snp_list} "
  1187. f"--max-maf {max_maf} --min-mac {min_mac} "
  1188. f"--wind {wind} --out {out}")
  1189. return gcta_run(args, work_dir)
  1192. # ============================================================================
  1193. # File management tools
  1194. # ============================================================================
  1196. @mcp.tool()
  1197. def gcta_list_files(directory: str = ".", pattern: str = "*") -> str:
  1198. """List files in a directory, optionally filtered by pattern.
  1200. Args:
  1201. directory: Directory path (default: current directory).
  1202. pattern: Glob pattern (default: "*" = all files).
  1203. Examples: "*.log", "*.hsq", "*.grm.*", "test*"
  1205. Returns:
  1206. JSON with list of files and their sizes.
  1207. """
  1208. if not os.path.exists(directory):
  1209. return json.dumps({"error": f"Directory not found: {directory}"}, indent=2)
  1211. files = []
  1212. search_pattern = os.path.join(directory, pattern)
  1213. for f in sorted(glob.glob(search_pattern)):
  1214. if os.path.isfile(f):
  1215. files.append({
  1216. "file": os.path.basename(f),
  1217. "path": os.path.abspath(f),
  1218. "size_bytes": os.path.getsize(f),
  1219. })
  1221. return json.dumps({
  1222. "directory": os.path.abspath(directory),
  1223. "pattern": pattern,
  1224. "file_count": len(files),
  1225. "files": files,
  1226. }, indent=2, ensure_ascii=False)
  1229. @mcp.tool()
  1230. def gcta_read_file(filepath: str, max_lines: int = 200) -> str:
  1231. """Read and return the content of a text file.
  1233. Useful for inspecting GCTA output files such as .log, .hsq, .mlma,
  1234. .ma, etc.
  1236. Args:
  1237. filepath: Path to the file to read.
  1238. max_lines: Maximum number of lines to return (default 200).
  1240. Returns:
  1241. JSON with file content (truncated to max_lines).
  1242. """
  1243. if not os.path.exists(filepath):
  1244. return json.dumps({"error": f"File not found: {filepath}"}, indent=2)
  1246. try:
  1247. with open(filepath, "r", errors="replace") as f:
  1248. lines = []
  1249. for i, line in enumerate(f):
  1250. if i >= max_lines:
  1251. lines.append(f"... (truncated at {max_lines} lines)")
  1252. break
  1253. lines.append(line.rstrip("\n\r"))
  1254. return json.dumps({
  1255. "file": os.path.abspath(filepath),
  1256. "lines_shown": len(lines),
  1257. "content": "\n".join(lines),
  1258. }, indent=2, ensure_ascii=False)
  1259. except Exception as e:
  1260. return json.dumps({"error": str(e)}, indent=2)
  1263. # ============================================================================
  1264. # Help documentation
  1265. # ============================================================================
  1267. def _get_help_docs() -> dict:
  1268. """Return a dictionary of help documentation by topic."""
  1269. return {
  1270. "all": _help_all(),
  1271. "data": _help_data(),
  1272. "grm": _help_grm(),
  1273. "reml": _help_reml(),
  1274. "bivariate_reml": _help_bivariate_reml(),
  1275. "hereg": _help_hereg(),
  1276. "pca": _help_pca(),
  1277. "mlma": _help_mlma(),
  1278. "cojo": _help_cojo(),
  1279. "gsmr": _help_gsmr(),
  1280. "mtcojo": _help_mtcojo(),
  1281. "fastgwa": _help_fastgwa(),
  1282. "fastbat": _help_fastbat(),
  1283. "simulation": _help_simulation(),
  1284. "ld": _help_ld(),
  1285. "fst": _help_fst(),
  1286. "expression": _help_expression(),
  1287. "acat": _help_acat(),
  1288. "options": _help_options(),
  1289. }
  1292. def _help_all():
  1293. return """GCTA (Genome-wide Complex Trait Analysis) - Available Analyses
  1295. GCTA is a tool for genome-wide association study (GWAS) analyses. It supports:
  1297. 1. DATA MANAGEMENT (gcta_make_bed, gcta_run)
  1298. - Read PLINK binary PED format (.bed/.bim/.fam)
  1299. - Filter SNPs and individuals
  1300. - Output filtered data in PLINK bed format
  1301. - Calculate allele frequencies, recode genotypes
  1303. 2. GRM CONSTRUCTION (gcta_make_grm)
  1304. - Calculate Genetic Relationship Matrix (GRM) from genotype data
  1305. - Support for additive, dominance, X-chromosome, and inbreeding GRMs
  1306. - Multi-component GRM analysis
  1308. 3. REML / GREML (gcta_reml)
  1309. - Estimate SNP-based heritability (h2_SNP) using REML
  1310. - Support for covariates, prevalence, and prior specification
  1311. - BLUP prediction of random effects
  1313. 4. BIVARIATE REML (gcta_bivariate_reml)
  1314. - Estimate genetic correlation (rg) between two traits
  1315. - Bivariate GREML analysis
  1317. 5. HE REGRESSION (gcta_hereg)
  1318. - Haseman-Elston regression for variance component estimation
  1319. - Faster but less precise than REML
  1321. 6. PCA (gcta_pca)
  1322. - Principal component analysis from GRM
  1323. - Population structure analysis
  1325. 7. MLMA (gcta_mlma)
  1326. - Mixed Linear Model Association analysis
  1327. - MLMA-LOCO (Leave-One-Chromosome-Out) for efficiency
  1328. - Genome-wide association testing
  1330. 8. COJO (gcta_cojo)
  1331. - Conditional and joint analysis of GWAS summary data
  1332. - Stepwise selection, conditional analysis, joint analysis
  1333. - SBLUP (summary-data BLUP)
  1335. 9. GSMR (gcta_gsmr)
  1336. - Generalized Summary-data Mendelian Randomization
  1337. - Causal effect estimation using GWAS summary data
  1339. 10. mtCOJO (gcta_mtcojo)
  1340. - Multi-trait COJO analysis
  1341. - Conditional analysis across multiple traits
  1343. 11. fastGWA (gcta_fastgwa)
  1344. - Efficient mixed-model association analysis
  1345. - fastGWA-MLM, fastGWA-MLM-exact, fastGWA-LR methods
  1346. - Sparse GRM for improved speed
  1348. 12. fastBAT (gcta_fastbat)
  1349. - Gene-based association test
  1350. - Combines SNP-level p-values within genes
  1352. 13. SIMULATION (gcta_simu_qt, gcta_simu_cc)
  1353. - Simulate quantitative trait phenotypes
  1354. - Simulate case-control phenotypes
  1355. - Based on real genotype data
  1357. 14. LD ANALYSIS (gcta_ld_pruning, gcta_ld_score)
  1358. - LD pruning to identify independent SNPs
  1359. - LD score calculation
  1361. 15. FST (gcta_fst)
  1362. - Fixation index for population differentiation analysis
  1364. 16. ACAT (gcta_acat)
  1365. - Aggregated Cauchy Association Test
  1366. - Gene-based test combining SNP p-values
  1368. 17. EXPRESSION DATA (gcta_run)
  1369. - Expression-based relationship matrix (ERM)
  1370. - ecojo analysis
  1371. - make-erm, make-erm-gz
  1373. Use gcta_help(topic) for detailed documentation on any analysis type.
  1374. Use gcta_run(args) for full command-line access to all GCTA options.
  1375. """
  1378. def _help_data():
  1379. return """GCTA Data Management
  1381. INPUT FORMAT:
  1382. --bfile <prefix> PLINK binary PED format (.bed/.bim/.fam)
  1383. --mbfile <file> List of multiple bfile prefixes
  1384. --bfile2 <prefix> Second dataset for comparison
  1385. --dosage-mach <dose> <info> Mach dosage format
  1386. --dosage-mach-gz <dose> <info> Mach dosage (gzipped)
  1387. --dosage-beagle <dose> <info> Beagle dosage format
  1389. FILTERING:
  1390. --maf <val> Minor allele frequency filter (0-0.5)
  1391. --max-maf <val> Maximum MAF filter
  1392. --chr <num> Restrict to chromosome
  1393. --autosome Restrict to autosomes
  1394. --autosome-num <num> Number of autosomes (default 22)
  1395. --extract <file> Include only SNPs in file
  1396. --exclude <file> Exclude SNPs in file
  1397. --keep <file> Keep only individuals in file
  1398. --remove <file> Remove individuals in file
  1400. OUTPUT:
  1401. --make-bed Output in PLINK bed format
  1402. --freq Calculate and output allele frequencies
  1403. --recode Recode genotypes
  1404. --recode-nomiss Recode without missing data
  1405. --recode-std Recode with standardization
  1406. --out <prefix> Output file prefix
  1408. EXAMPLES:
  1409. # Filter SNPs by MAF and output PLINK bed
  1410. gcta_run("--bfile test --make-bed --maf 0.05 --out filtered")
  1412. # Extract chromosome 1 only
  1413. gcta_run("--bfile test --make-bed --chr 1 --out chr1")
  1415. # Calculate allele frequencies
  1416. gcta_run("--bfile test --freq --out freq")
  1417. """
  1420. def _help_grm():
  1421. return """GCTA GRM (Genetic Relationship Matrix) Construction
  1423. The GRM estimates pairwise genomic relatedness between individuals based
  1424. on genome-wide SNPs. It is the foundational step for most downstream
  1425. analyses (REML, PCA, MLMA).
  1427. BASIC COMMAND:
  1428. --bfile <prefix> --make-grm --out <output_prefix>
  1430. OUTPUT FILES:
  1431. {prefix}.grm.bin - GRM binary data (lower triangle, column-major)
  1432. {prefix}.grm.N.bin - Number of SNPs used (same format)
  1433. {prefix}.grm.id - Individual IDs (FID, IID)
  1435. OPTIONS:
  1436. --make-grm-alg <0|1> Algorithm: 0 = VanRaden (default), 1 = alternative
  1437. --grm-adj <val> GRM adjustment factor (0-1)
  1438. --grm-cutoff <val> Remove pairs with GRM > cutoff
  1439. --dc <0|1> Dosage compensation (0=none, 1=X chr)
  1440. --dominance Calculate dominance GRM
  1441. --make-grm-xchr Calculate X-chromosome GRM
  1442. --make-grm-inbred Calculate inbreeding GRM
  1443. --save-ram Use memory-saving mode
  1444. --autosome-num <num> Number of autosomes
  1445. --autosome Restrict to autosomes
  1446. --chr <num> Restrict to specific chromosome
  1447. --maf <val> MAF filter
  1448. --thread-num <num> Number of CPU threads
  1450. MULTI-GRM:
  1451. --mgrm <file> Multiple GRM file list
  1452. --mgrm-bin <file> Multiple GRM binary file list
  1453. --mgrm-gz <file> Multiple GRM gzipped file list
  1455. EXAMPLES:
  1456. # Basic GRM
  1457. gcta_make_grm(bfile="test", out="test_grm")
  1459. # GRM with MAF filter and 4 threads
  1460. gcta_make_grm(bfile="test", out="test_grm", maf=0.01, thread_num=4)
  1462. # Dominance GRM
  1463. gcta_make_grm(bfile="test", out="test_dgrm", dominance=True)
  1464. """
  1467. def _help_reml():
  1468. return """GCTA REML / GREML Analysis
  1470. REML (Restricted Maximum Likelihood) estimates the proportion of
  1471. phenotypic variance explained by all genome-wide SNPs (SNP-based
  1472. heritability, h2_SNP).
  1474. BASIC COMMAND:
  1475. --grm <prefix> --pheno <file> --reml --out <output_prefix>
  1477. REQUIRED INPUT:
  1478. --grm <prefix> GRM file prefix (from gcta_make_grm)
  1479. --pheno <file> Phenotype file (FID, IID, phenotype)
  1480. --reml Perform REML analysis
  1481. --out <prefix> Output file prefix
  1483. OUTPUT FILES:
  1484. {prefix}.hsq - Heritability estimates and SE
  1485. {prefix}.log - Analysis log
  1487. COVARIATES:
  1488. --qcovar <file> Quantitative covariate file
  1489. --covar <file> Categorical covariate file
  1491. REML OPTIONS:
  1492. --reml-maxit <num> Maximum iterations (default 100)
  1493. --reml-priors <vals> Prior values of variance explained
  1494. --reml-priors-var <vals> Prior values of variance components
  1495. --reml-no-constrain Do not constrain variance components
  1496. --reml-no-lrt Skip likelihood ratio test
  1497. --reml-pred-rand Predict random effects (BLUP)
  1498. --reml-est-fix Estimate fixed effects
  1500. OTHER OPTIONS:
  1501. --prevalence <val> Disease prevalence for liability conversion
  1502. --grm-cutoff <val> GRM cutoff for related pairs
  1503. --gxe <file> GxE interaction file
  1504. --mpheno <num> Phenotype column number
  1505. --thread-num <num> Number of CPU threads
  1507. EXAMPLES:
  1508. # Basic REML
  1509. gcta_reml(grm="test_grm", pheno="test.phen", out="test_reml")
  1511. # REML with covariates
  1512. gcta_reml(grm="test_grm", pheno="test.phen", out="test_reml",
  1513. qcovar="covariates.txt", covar="sex.txt")
  1514. """
  1517. def _help_bivariate_reml():
  1518. return """GCTA Bivariate REML Analysis
  1520. Bivariate REML estimates the genetic correlation (rg) between two traits
  1521. using bivariate GREML.
  1523. BASIC COMMAND:
  1524. --grm <prefix> --pheno <file> --reml-bivar <mphen> <mphen2> --out <prefix>
  1526. REQUIRED INPUT:
  1527. --grm <prefix> GRM file prefix
  1528. --pheno <file> Phenotype file (with 2+ trait columns)
  1529. --reml-bivar <p1> <p2> Bivariate REML with trait columns p1 and p2
  1530. --out <prefix> Output file prefix
  1532. OUTPUT FILES:
  1533. {prefix}.hsq - Variance components and genetic correlation
  1535. OPTIONS:
  1536. --reml-maxit <num> Maximum iterations
  1537. --reml-bivar-prevalence <k1> <k2> Prevalence for both traits
  1538. --reml-bivar-nocove Ignore residual covariance
  1539. --reml-bivar-no-constrain No constrain on variance components
  1540. --qcovar <file> Quantitative covariate file
  1541. --covar <file> Categorical covariate file
  1542. --grm-cutoff <val> GRM cutoff
  1543. --thread-num <num> Number of CPU threads
  1545. EXAMPLE:
  1546. gcta_bivariate_reml(grm="test_grm", pheno="traits.phen",
  1547. out="bireml", mphen=1, mphen2=2)
  1548. """
  1551. def _help_hereg():
  1552. return """GCTA Haseman-Elston (HE) Regression
  1554. HE regression is a method-of-moments estimator for variance components.
  1555. It is faster but less precise than REML. Useful as a sanity check or
  1556. when REML fails to converge.
  1558. BASIC COMMAND:
  1559. --grm <prefix> --pheno <file> --HEreg --out <prefix>
  1561. OUTPUT FILES:
  1562. {prefix}.hsq - HE regression estimates
  1564. OPTIONS:
  1565. --qcovar <file> Quantitative covariate file
  1566. --covar <file> Categorical covariate file
  1567. --mpheno <num> Phenotype column number
  1568. --grm-cutoff <val> GRM cutoff
  1569. --thread-num <num> Number of CPU threads
  1571. EXAMPLE:
  1572. gcta_hereg(grm="test_grm", pheno="test.phen", out="test_he")
  1573. """
  1576. def _help_pca():
  1577. return """GCTA Principal Component Analysis (PCA)
  1579. PCA from the GRM captures population structure. Output eigenvalues
  1580. and eigenvectors are saved.
  1582. BASIC COMMAND:
  1583. --grm <prefix> --pca <num_pcs> --out <prefix>
  1585. OUTPUT FILES:
  1586. {prefix}.eigenval - Eigenvalues (one per line)
  1587. {prefix}.eigenvec - Eigenvectors (FID, IID, PC1, PC2, ...)
  1589. OPTIONS:
  1590. --pca <num> Number of PCs (default 20)
  1591. --thread-num <num> Number of CPU threads
  1593. EXAMPLE:
  1594. gcta_pca(grm="test_grm", out="test_pca", pc_num=10)
  1595. """
  1598. def _help_mlma():
  1599. return """GCTA Mixed Linear Model Association (MLMA)
  1601. MLMA performs genome-wide association analysis using a linear mixed
  1602. model that accounts for population structure via the GRM.
  1604. METHODS:
  1605. --mlma Standard MLMA (uses full GRM)
  1606. --mlma-loco MLMA with Leave-One-Chromosome-Out (faster, recommended)
  1608. BASIC COMMAND:
  1609. --bfile <prefix> --grm <prefix> --pheno <file> --mlma-loco --out <prefix>
  1611. REQUIRED INPUT:
  1612. --bfile <prefix> PLINK binary file prefix (genotype data)
  1613. --grm <prefix> GRM file prefix
  1614. --pheno <file> Phenotype file
  1615. --out <prefix> Output file prefix
  1617. OUTPUT FILES:
  1618. {prefix}.mlma - Association results (Chr, SNP, BP, A1, F, BETA, SE, P)
  1620. COVARIATES:
  1621. --qcovar <file> Quantitative covariate file
  1622. --covar <file> Categorical covariate file
  1624. OPTIONS:
  1625. --mpheno <num> Phenotype column number
  1626. --reml-maxit <num> Maximum REML iterations
  1627. --maf <val> Minor allele frequency filter
  1628. --autosome-num <num> Number of autosomes
  1629. --thread-num <num> Number of CPU threads
  1631. EXAMPLES:
  1632. # MLMA-LOCO (recommended)
  1633. gcta_mlma(bfile="test", grm="test_grm", pheno="test.phen",
  1634. out="test_mlma", loco=True)
  1636. # Standard MLMA
  1637. gcta_mlma(bfile="test", grm="test_grm", pheno="test.phen",
  1638. out="test_mlma", loco=False)
  1639. """
  1642. def _help_cojo():
  1643. return """GCTA Conditional and Joint (COJO) Analysis
  1645. COJO performs stepwise selection, conditional analysis, or joint analysis
  1646. of GWAS summary statistics, using an LD reference panel from PLINK bed
  1647. data.
  1649. BASIC COMMAND:
  1650. --bfile <prefix> --cojo-file <file> --cojo-slct --out <prefix>
  1652. REQUIRED INPUT:
  1653. --bfile <prefix> PLINK binary file prefix (LD reference panel)
  1654. --cojo-file <file> GWAS summary data file (.ma format)
  1655. --out <prefix> Output file prefix
  1657. GWAS SUMMARY FILE FORMAT (.ma):
  1658. Columns: SNP A1 A2 freq b se p n
  1659. (Header line required)
  1661. METHODS:
  1662. --cojo-slct Stepwise selection (backward + forward)
  1663. --cojo-forward Forward-only selection
  1664. --cojo-backward Backward-only selection
  1665. --cojo-cond <file> Conditional analysis (condition on given SNPs)
  1666. --cojo-joint Joint analysis (all SNPs in one model)
  1667. --cojo-sblup <factor> SBLUP (summary-data BLUP)
  1669. OPTIONS:
  1670. --cojo-p <val> P-value threshold (default 5e-8)
  1671. --cojo-wind <kb> Window size in Kb (default 10000)
  1672. --cojo-collinear <val> Collinearity threshold (default 0.9)
  1673. --cojo-gc <val> Genomic control inflation factor
  1674. --cojo-top-SNPs <num> Number of top SNPs to select
  1675. --maf <val> Minor allele frequency filter
  1676. --thread-num <num> Number of CPU threads
  1678. OUTPUT FILES:
  1679. {prefix}.jma.cojo - Joint analysis results
  1680. {prefix}.cma.cojo - Conditional analysis results
  1681. {prefix}.log - Analysis log
  1683. EXAMPLES:
  1684. # Stepwise selection
  1685. gcta_cojo(bfile="test", cojo_file="gwas.ma", out="cojo_res",
  1686. method="slct")
  1688. # Conditional analysis
  1689. gcta_cojo(bfile="test", cojo_file="gwas.ma", out="cojo_res",
  1690. method="cond", cojo_cond_snplist="cond_snps.txt")
  1691. """
  1694. def _help_gsmr():
  1695. return """GCTA GSMR (Generalized Summary-data Mendelian Randomization)
  1697. GSMR estimates the causal effect of an exposure on an outcome using
  1698. GWAS summary data and an LD reference panel.
  1700. BASIC COMMAND:
  1701. --bfile <prefix> --gsmr-file <file> --gsmr-direction <0|1|2> --out <prefix>
  1703. REQUIRED INPUT:
  1704. --bfile <prefix> PLINK binary file prefix (LD reference)
  1705. --gsmr-file <file> File listing exposure and outcome GWAS files
  1706. --gsmr-direction <num> 0 = forward, 1 = reverse, 2 = bi-directional
  1707. --out <prefix> Output file prefix
  1709. GSMR FILE FORMAT:
  1710. Exposure_GWAS_filename Outcome_GWAS_filename n_exp n_out
  1712. OPTIONS:
  1713. --gsmr2-beta Use GSMR2-beta version (multi-SNP HEIDI)
  1714. --gsmr-snp-min <num> Minimum number of SNP instruments
  1715. --gwas-thresh <val> GWAS p-value threshold
  1716. --clump-kb <val> Clumping window size
  1717. --clump-r2 <val> Clumping r2 threshold
  1718. --heidi-thresh <val> HEIDI-outlier p-value threshold
  1719. --diff-freq <val> Frequency difference threshold
  1720. --thread-num <num> Number of CPU threads
  1722. OUTPUT FILES:
  1723. {prefix}.gsmr - GSMR results
  1724. {prefix}.log - Analysis log
  1726. EXAMPLE:
  1727. gcta_gsmr(bfile="test", gsmr_file="gsmr_data.txt",
  1728. out="gsmr_res", gsmr_direction=0)
  1729. """
  1732. def _help_mtcojo():
  1733. return """GCTA mtCOJO (Multi-trait COJO) Analysis
  1735. mtCOJO performs conditional analysis across multiple traits using GWAS
  1736. summary data, adjusting for the genetic correlation between traits.
  1738. BASIC COMMAND:
  1739. --bfile <prefix> --mtcojo-file <file> --mtcojo-bxy <file> --out <prefix>
  1741. REQUIRED INPUT:
  1742. --bfile <prefix> PLINK binary file prefix (LD reference)
  1743. --mtcojo-file <file> File listing trait GWAS summary files
  1744. --mtcojo-bxy <file> File with causal effect estimates between traits
  1745. --out <prefix> Output file prefix
  1747. OPTIONS:
  1748. --gwas-thresh <val> GWAS p-value threshold
  1749. --clump-kb <val> Clumping window size
  1750. --clump-r2 <val> Clumping r2 threshold
  1751. --diff-freq <val> Frequency difference threshold
  1752. --heidi-thresh <val> HEIDI-outlier threshold
  1753. --thread-num <num> Number of CPU threads
  1755. EXAMPLE:
  1756. gcta_mtcojo(bfile="test", mtcojo_file="traits.txt",
  1757. mtcojo_bxy="bxy.txt", out="mtcojo_res")
  1758. """
  1761. def _help_fastgwa():
  1762. return """GCTA fastGWA Analysis
  1764. fastGWA is an efficient mixed-model association method that uses sparse
  1765. GRM techniques for improved speed on large datasets.
  1767. METHODS:
  1768. --fastGWA fastGWA (requires --grm sparse GRM)
  1769. --fastGWA-mlm fastGWA-MLM (linear mixed model, recommended)
  1770. --fastGWA-mlm-exact fastGWA-MLM exact mode
  1771. --fastGWA-lr fastGWA linear regression (no GRM)
  1773. BASIC COMMAND:
  1774. --bfile <prefix> --pheno <file> --fastGWA-mlm --out <prefix>
  1776. REQUIRED INPUT:
  1777. --bfile <prefix> PLINK binary file prefix
  1778. --pheno <file> Phenotype file
  1779. --out <prefix> Output file prefix
  1781. COVARIATES:
  1782. --qcovar <file> Quantitative covariate file
  1783. --covar <file> Categorical covariate file
  1785. OPTIONS:
  1786. --grm <prefix> GRM file prefix (for fastGWA with GRM)
  1787. --maf <val> Minor allele frequency filter
  1788. --autosome-num <num> Number of autosomes
  1789. --thread-num <num> Number of CPU threads
  1791. OUTPUT FILES:
  1792. {prefix}.fastGWA - Association results
  1793. {prefix}.log - Analysis log
  1795. EXAMPLES:
  1796. # fastGWA-MLM (recommended)
  1797. gcta_fastgwa(bfile="test", pheno="test.phen", out="fastgwa_res",
  1798. method="fastGWA-mlm")
  1800. # fastGWA-LR (no GRM, fastest)
  1801. gcta_fastgwa(bfile="test", pheno="test.phen", out="fastgwa_res",
  1802. method="fastGWA-lr")
  1803. """
  1806. def _help_fastbat():
  1807. return """GCTA fastBAT Gene-based Association Test
  1809. fastBAT performs gene-based association tests by combining SNP-level
  1810. p-values within genes, accounting for LD structure.
  1812. BASIC COMMAND:
  1813. --bfile <prefix> --fastBAT <gwas_summary> --out <prefix>
  1815. REQUIRED INPUT:
  1816. --bfile <prefix> PLINK binary file prefix (LD reference)
  1817. --fastBAT <file> GWAS summary data file (.ma format)
  1818. --fastBAT-gene-list <file> Gene annotation file
  1819. --fastBAT-set-list <file> SNP set file (predefined sets)
  1820. --out <prefix> Output file prefix
  1822. GWAS SUMMARY FILE FORMAT (.ma):
  1823. Columns: SNP A1 A2 freq b se p n
  1825. OPTIONS:
  1826. --fastBAT-wind <kb> Window size in Kb (default 50)
  1827. --fastBAT-ld-cutoff <val> LD r2 cutoff (default 0.9)
  1828. --fastBAT-write-snpset Write SNP set file
  1829. --maf <val> Minor allele frequency filter
  1830. --thread-num <num> Number of CPU threads
  1832. OUTPUT FILES:
  1833. {prefix}.fastBAT - Gene-based test results
  1834. {prefix}.log - Analysis log
  1836. EXAMPLE:
  1837. gcta_fastbat(bfile="test", fastbat_file="gwas.ma",
  1838. out="fastbat_res", fastbat_gene_list="genes.txt")
  1839. """
  1842. def _help_simulation():
  1843. return """GCTA Phenotype Simulation
  1845. GCTA can simulate phenotypes based on real genotype data, useful for
  1846. power analysis and method evaluation.
  1848. QUANTITATIVE TRAIT SIMULATION:
  1849. --bfile <prefix> --simu-qt --simu-hsq <h2> --out <prefix>
  1851. CASE-CONTROL SIMULATION:
  1852. --bfile <prefix> --simu-cc <cases> <controls> --simu-hsq <h2>
  1853. --simu-k <prevalence> --out <prefix>
  1855. OPTIONS:
  1856. --simu-hsq <val> Simulated heritability (default 0.1)
  1857. --simu-k <val> Disease prevalence (default 0.1)
  1858. --simu-rep <num> Number of repetitions (default 1)
  1859. --simu-causal-loci <file> File of causal loci
  1860. --simu-seed <val> Random seed
  1861. --simu-eff-mod <0|1> Effect model: 0 = additive, 1 = non-additive
  1862. --thread-num <num> Number of CPU threads
  1864. EXAMPLES:
  1865. # Simulate QT with h2=0.5
  1866. gcta_simu_qt(bfile="test", out="simu_qt", simu_hsq=0.5)
  1868. # Simulate CC with 1000 cases, 1000 controls
  1869. gcta_simu_cc(bfile="test", simu_case_num=1000, simu_control_num=1000,
  1870. out="simu_cc", simu_hsq=0.3, simu_k=0.05)
  1871. """
  1874. def _help_ld():
  1875. return """GCTA LD Analysis
  1877. GCTA provides several LD analysis tools:
  1879. LD PRUNING:
  1880. --bfile <prefix> --ld-pruning <r2> --ld-wind <kb> --out <prefix>
  1881. Removes SNPs in high LD, retaining approximately independent SNPs.
  1883. LD SCORE CALCULATION:
  1884. --bfile <prefix> --ld-score --ld-wind <kb> --out <prefix>
  1885. Computes LD scores (sum of LD r2 with neighboring SNPs).
  1887. LD SCORE REGRESSION:
  1888. --bfile <prefix> --ld-score-region --ld-wind <kb> --out <prefix>
  1889. Performs LD score regression to estimate heritability and confounding.
  1891. OPTIONS:
  1892. --ld-pruning <r2> LD r2 threshold for pruning
  1893. --ld-wind <kb> LD window size in Kb
  1894. --ld-score Calculate LD scores
  1895. --ld-score-region LD score regression
  1896. --ld-score-multi <f> Multi-set LD score
  1897. --ld-rsq-cutoff <val> LD r2 cutoff
  1898. --ld-max-rsq Maximum LD r2
  1899. --ld-step <num> LD step size
  1900. --thread-num <num> Number of CPU threads
  1902. EXAMPLES:
  1903. # LD pruning with r2=0.1
  1904. gcta_ld_pruning(bfile="test", ld_pruning_rsq=0.1, out="pruned")
  1906. # LD score calculation
  1907. gcta_ld_score(bfile="test", out="ldscores")
  1908. """
  1911. def _help_fst():
  1912. return """GCTA Fst (Fixation Index) Analysis
  1914. Fst measures population differentiation. GCTA computes pairwise Fst
  1915. between subpopulations.
  1917. BASIC COMMAND:
  1918. --bfile <prefix> --fst --sub-popu <file> --out <prefix>
  1920. REQUIRED INPUT:
  1921. --bfile <prefix> PLINK binary file prefix
  1922. --sub-popu <file> Subpopulation file (FID, IID, pop ID)
  1923. --out <prefix> Output file prefix
  1925. OPTIONS:
  1926. --maf <val> Minor allele frequency filter
  1927. --thread-num <num> Number of CPU threads
  1929. OUTPUT FILES:
  1930. {prefix}.Fst - Fst results
  1931. {prefix}.log - Analysis log
  1933. EXAMPLE:
  1934. gcta_fst(bfile="test", sub_popu="pops.txt", out="fst_res")
  1935. """
  1938. def _help_expression():
  1939. return """GCTA Expression Data Analysis
  1941. GCTA supports analysis of gene expression data:
  1943. EXPRESSION FILE INPUT:
  1944. --efile <file> Expression data file
  1946. EXPRESSION RELATIONSHIP MATRIX (ERM):
  1947. --make-erm Make ERM (binary output)
  1948. --make-erm-gz Make ERM (gzipped text output)
  1949. --make-erm-alg <1|2|3> Algorithm for ERM calculation
  1951. ECOJO ANALYSIS:
  1952. --ecojo <file> ecojo analysis with MA file
  1953. --ecojo-slct ecojo stepwise selection
  1954. --ecojo-p <val> ecojo p-value threshold
  1955. --ecojo-collinear <v> ecojo collinearity threshold
  1956. --ecojo-blup <lambda> ecojo BLUP with lambda
  1958. E-COR (Expression Correlation):
  1959. --e-cor <file> Expression correlation file
  1961. EXAMPLE:
  1962. gcta_run("--efile expr.txt --make-erm --out expr_erm")
  1963. """
  1966. def _help_acat():
  1967. return """GCTA ACAT (Aggregated Cauchy Association Test)
  1969. ACAT combines p-values from multiple SNPs within a gene using a Cauchy
  1970. combination method, providing a gene-level association test.
  1972. BASIC COMMAND:
  1973. --acat --gene-list <file> --snp-list <file> --out <file>
  1975. REQUIRED INPUT:
  1976. --gene-list <file> Gene list file (gene, chr, start, end)
  1977. --snp-list <file> SNP list file with p-values
  1978. --out <file> Output file (default acat_res.csv)
  1980. OPTIONS:
  1981. --max-maf <val> Maximum MAF for included SNPs (default 0.01)
  1982. --min-mac <num> Minimum minor allele count (default 20)
  1983. --wind <num> Extension length for gene boundaries
  1985. OUTPUT FILE:
  1986. {out} - ACAT gene-based test results
  1988. EXAMPLE:
  1989. gcta_acat(gene_list="genes.txt", snp_list="snps.txt",
  1990. out="acat_results.csv")
  1991. """
  1994. def _help_options():
  1995. return """GCTA Complete Command-Line Options Reference
  1997. DATA INPUT:
  1998. --bfile <prefix> PLINK binary PED format
  1999. --mbfile <file> Multiple bfile list
  2000. --bfile2 <prefix> Second bfile
  2001. --dosage-mach <d> <i> Mach dosage
  2002. --dosage-mach-gz <d> <i> Mach dosage (gzipped)
  2003. --dosage-beagle <d> <i> Beagle dosage
  2005. DATA MANAGEMENT:
  2006. --make-bed Output PLINK bed format
  2007. --freq / --freqx Calculate allele frequencies
  2008. --recode / --recode-nomiss / --recode-std Recode genotypes
  2009. --keep <file> Keep individuals
  2010. --remove <file> Remove individuals
  2011. --extract <file> Include SNPs
  2012. --exclude <file> Exclude SNPs
  2013. --chr <num> Chromosome filter
  2014. --autosome Autosome filter
  2015. --autosome-num <num> Number of autosomes
  2016. --maf <val> MAF filter
  2017. --max-maf <val> Max MAF filter
  2018. --update-sex <file> Update sex
  2019. --update-freq <file> Update allele frequencies
  2020. --update-ref-allele <f> Update reference allele
  2021. --save-ram Memory-saving mode
  2023. GRM:
  2024. --make-grm Make GRM (binary output)
  2025. --make-grm-gz Make GRM (text output)
  2026. --make-grm-alg <0|1> GRM algorithm
  2027. --make-grm-d Dominance GRM
  2028. --make-grm-xchr X-chromosome GRM
  2029. --make-grm-inbred Inbreeding GRM
  2030. --grm <prefix> Read GRM
  2031. --grm-gz <prefix> Read GRM (gzipped)
  2032. --mgrm <file> Multiple GRMs
  2033. --mgrm-gz <file> Multiple GRMs (gzipped)
  2034. --grm-cutoff <val> GRM cutoff
  2035. --grm-adj <val> GRM adjustment
  2036. --dc <0|1> Dosage compensation
  2038. PCA:
  2039. --pca <num> PCA with num PCs
  2040. --pc-loading <file> PC loading
  2041. --project-loading <f> <n> Project loading
  2043. REML:
  2044. --reml REML analysis
  2045. --reml-bivar <p1> <p2> Bivariate REML
  2046. --reml-maxit <num> Max iterations
  2047. --reml-priors <vals> Prior variance explained
  2048. --reml-priors-var <vals> Prior variance components
  2049. --reml-no-constrain No constrain
  2050. --reml-no-lrt No LRT
  2051. --reml-pred-rand Predict random effects
  2052. --reml-est-fix Estimate fixed effects
  2053. --prevalence <val> Disease prevalence
  2055. HE REGRESSION:
  2056. --HEreg HE regression
  2057. --HEreg-bivar <p1> <p2> Bivariate HE regression
  2059. MLMA:
  2060. --mlma MLMA association
  2061. --mlma-loco MLMA-LOCO
  2062. --mlma-no-adj-covar No adj covar
  2064. COJO:
  2065. --cojo-file <file> GWAS summary data
  2066. --cojo-slct Stepwise selection
  2067. --cojo-forward Forward selection
  2068. --cojo-backward Backward selection
  2069. --cojo-joint Joint analysis
  2070. --cojo-cond <file> Conditional analysis
  2071. --cojo-sblup <factor> SBLUP
  2072. --cojo-p <val> P-value threshold
  2073. --cojo-wind <kb> Window size
  2074. --cojo-collinear <val> Collinearity threshold
  2075. --cojo-gc <val> Genomic control
  2076. --cojo-top-SNPs <num> Top SNPs
  2078. GSMR:
  2079. --gsmr-file <file> GSMR data file
  2080. --gsmr-direction <0|1|2> Direction
  2081. --gsmr2-beta GSMR2-beta
  2082. --gsmr-snp-min <num> Min SNP instruments
  2083. --gwas-thresh <val> GWAS threshold
  2084. --heidi-thresh <val> HEIDI threshold
  2086. mtCOJO:
  2087. --mtcojo-file <file> mtCOJO data file
  2088. --mtcojo-bxy <file> Bxy file
  2090. fastGWA:
  2091. --fastGWA fastGWA
  2092. --fastGWA-mlm fastGWA-MLM
  2093. --fastGWA-mlm-exact fastGWA-MLM exact
  2094. --fastGWA-lr fastGWA-LR
  2096. fastBAT:
  2097. --fastBAT <file> fastBAT GWAS summary
  2098. --fastBAT-gene-list <f> Gene list
  2099. --fastBAT-set-list <f> SNP set list
  2100. --fastBAT-wind <kb> Window size
  2101. --fastBAT-ld-cutoff <v> LD cutoff
  2103. SIMULATION:
  2104. --simu-qt Simulate QT
  2105. --simu-cc <cases> <ctrl> Simulate CC
  2106. --simu-hsq <val> Heritability
  2107. --simu-k <val> Prevalence
  2108. --simu-rep <num> Repetitions
  2109. --simu-causal-loci <file> Causal loci
  2110. --simu-seed <val> Seed
  2112. LD:
  2113. --ld <file> LD analysis
  2114. --ld-pruning <r2> LD pruning
  2115. --ld-score LD score
  2116. --ld-score-region LD score regression
  2117. --ld-wind <kb> LD window
  2118. --ld-rsq-cutoff <val> LD r2 cutoff
  2120. FST:
  2121. --fst Fst analysis
  2122. --sub-popu <file> Subpopulation file
  2124. EXPRESSION:
  2125. --efile <file> Expression file
  2126. --e-cor <file> Expression correlation
  2127. --make-erm Make ERM
  2128. --ecojo <file> ecojo analysis
  2130. ACAT:
  2131. --acat ACAT analysis
  2132. --gene-list <file> Gene list
  2133. --snp-list <file> SNP list
  2134. --max-maf <val> Max MAF
  2135. --min-mac <num> Min MAC
  2137. THREADING:
  2138. --thread-num <num> Number of threads
  2139. --threads <num> Alias for --thread-num
  2141. OUTPUT:
  2142. --out <prefix> Output prefix
  2143. """
  2146. # ============================================================================
  2147. # Main entry point
  2148. # ============================================================================
  2150. if __name__ == "__main__":
  2151. print(f"GCTA MCP Server starting...", file=sys.stderr)
  2152. print(f" Binary path: {GCTA_BINARY_PATH}", file=sys.stderr)
  2153. print(f" Work directory: {GCTA_WORK_DIR}", file=sys.stderr)
  2154. print(f" Timeout: {GCTA_TIMEOUT}s", file=sys.stderr)
  2155. mcp.run()

server.py at commit 2ff095c, under GPL-3.0 · at the source

Overview

Authors: Miguel Novoa-Bravo1, Jennifer RS Meadows2,3, Filipe Serra-Bragança4, Britt van de Vall5, Klas Kullander6, Marie Rhodin5, Gabriella Lindgren5
ORCID iDs: Miguel Novoa-Bravo, Britt van de Vall
  1. Genética Animal de Colombia SAS, Bogotá, Colombia
  2. Department of Medical Biochemistry and Microbiology, Uppsala University, 75123 Uppsala, Sweden
  3. SciLifeLab, Uppsala University, 75123 Uppsala, Sweden
  4. Department of Clinical Sciences, Faculty of Veterinary Medicine, Utrecht University, Utrecht 3584CM, the Netherlands
  5. Department of Animal Biosciences, Swedish University of Agricultural Sciences, P.O. Box 7023, 75007 Uppsala, Sweden
  6. Department of Immunology, Genetics and Pathology, Uppsala University, Uppsala, Sweden
Institutions: Uppsala University (Sweden); Science for Life Laboratory (Sweden); Utrecht University (Netherlands); Swedish University of Agricultural Sciences (Sweden)
Journal: iScience, volume 29, issue 7, article 116320
Dates: received 17 December 2025; accepted 25 May 2026; published online 11 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1016/j.isci.2026.116320 · PMID 42325571 · PMCID PMC13276308 · OpenAlex W7164341170
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), other (organism), developmental (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, fMRI & imaging, Physiology & signal measures
Keywords: genetics, developmental neuroscience, machine learning
Topic: Genetic Mapping and Diversity in Plants and Animals (Genetics, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Hjärnfonden (FO2022-0018); Vetenskapsrådet (2022-01245); Svenska Forskningsrådet Formas; Formas
Citations: not cited yet (Europe PMC); 86 references in the paper

Abstract

Coordinated mammalian locomotion relies on spinal circuits where DMRT3 regulates strides and alternative gaits. However, DMRT3 does not explain differences between the Colombian paso horse breed’s specialized gaits, Colombian trocha and Colombian trot. We used inertial sensors and machine learning to develop accurate phenotyping (n = 225 horses), before performing genome-wide association analysis (n = 85 horses, 670K array). We identified a 2.43 Mb quantitative trait locus (QTL) on ECA16, where haplotypes featuring lead variants (rs1147402472, p = 1.95 × 10−8; rs1136628503, p = 8.52 × 10−8) explained 48.6% of gait variance (p < 0.001). The QTL contained 11 genes linked to neurodevelopment and muscle regulation, including LHFPL4, SRGAP3, and ATP2B2. Interaction analyses suggested a functional link between the latter genes. rs1147402472-C was common across diverse gaited horse breeds, but the Colombian trot-specific haplotype (rs1147402472-C/rs1136628503-T) was absent elsewhere. These findings demonstrate that fine-scale genetic differentiation at ECA16 underlies neural adaptations distinguishing complex locomotor traits and highlight the power of AI-assisted phenotyping in genomics.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 4 matches between paragraphs and lines of code.

JianYang-Lab/GCTA

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 2ff095cad2ea5ca213c1b008401382aa7b05715e, 21 July 2026
Languages: C++ (43), C/C++ (26), Shell (3), Python (1), C (1)
Size: 106 files, 74 scripts
Software Heritage: not archived
Found in: the resources table
Holds: README, license file, environment (.opencode/gcta-mcp/pyproject.toml, .opencode/gcta-mcp/requirements.txt), documentation
Not found: CITATION.cff, tests, continuous integration
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
76 files
At the source: github.com/JianYang-Lab/GCTA

saigegit/SAIGE

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 365deeff9149fc2b5a66da8fa3799637bcb30af3, 25 June 2026
Languages: C++ (26), R (24), Shell (4), C/C++ (2)
Size: 481 files, 56 scripts
Software Heritage: not archived
Found in: the resources table
Holds: README, license file, environment (DESCRIPTION, pixi.toml, docker/Dockerfile, src/SAIGE/DESCRIPTION), documentation
Not found: CITATION.cff, tests, continuous integration
Tools: data.table (16 files), tidyverse (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
58 files
At the source: github.com/saigegit/SAIGE

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 130 scripts, each with its path and the digest of its content;
  • 4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data and code availability

Genotyping array data and metadata have been deposited at European Variation Archive (EVA): PRJEB105194, and are publicly available as of the date of publication.

Phenotype data have been deposited in Mendeley Dataset: https://doi.org/10.17632/52rzdhx297.1

This paper does not report original code.

Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

  • Authors: added Miguel Novoa-Bravo (0000-0002-9950-2590); removed Miguel Novoa-Bravo

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 3 keywords, 4 funders, 86 references.

Cite

This paper

Novoa-Bravo, M., Meadows, J. R., Serra-Bragança, F., van de Vall, B., Kullander, K., Rhodin, M., & Lindgren, G. (2026). Machine learning-based gait classification and genome-wide association identify a QTL for gait type in Colombian paso horses. iScience, 29(7), 116320. https://doi.org/10.1016/j.isci.2026.116320

BibTeX

@article{novoabravo2026machine,
author = {Novoa-Bravo, Miguel and Meadows, Jennifer RS and Serra-Bragança, Filipe and van de Vall, Britt and Kullander, Klas and Rhodin, Marie and Lindgren, Gabriella},
title = {{Machine learning-based gait classification and genome-wide association identify a QTL for gait type in Colombian paso horses}},
journal = {iScience},
year = {2026},
month = jun,
volume = {29},
number = {7},
pages = {116320},
publisher = {Elsevier},
issn = {2589-0042},
doi = {10.1016/j.isci.2026.116320},
url = {https://doi.org/10.1016/j.isci.2026.116320},
pmid = {42325571},
pmcid = {PMC13276308}
}

RIS

TY - JOUR
AU - Novoa-Bravo, Miguel
AU - Meadows, Jennifer RS
AU - Serra-Bragança, Filipe
AU - van de Vall, Britt
AU - Kullander, Klas
AU - Rhodin, Marie
AU - Lindgren, Gabriella
TI - Machine learning-based gait classification and genome-wide association identify a QTL for gait type in Colombian paso horses
T2 - iScience
J2 - iScience
PY - 2026
DA - 2026/06/11
VL - 29
IS - 7
SP - 116320
SN - 2589-0042
PB - Elsevier
DO - 10.1016/j.isci.2026.116320
UR - https://doi.org/10.1016/j.isci.2026.116320
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.isci.2026.116320",
"type": "article-journal",
"title": "Machine learning-based gait classification and genome-wide association identify a QTL for gait type in Colombian paso horses",
"container-title": "iScience",
"author": [
{
"family": "Novoa-Bravo",
"given": "Miguel"
},
{
"family": "Meadows",
"given": "Jennifer RS"
},
{
"family": "Serra-Bragança",
"given": "Filipe"
},
{
"family": "van de Vall",
"given": "Britt"
},
{
"family": "Kullander",
"given": "Klas"
},
{
"family": "Rhodin",
"given": "Marie"
},
{
"family": "Lindgren",
"given": "Gabriella"
}
],
"container-title-short": "iScience",
"volume": "29",
"issue": "7",
"page": "116320",
"DOI": "10.1016/j.isci.2026.116320",
"PMID": "42325571",
"PMCID": "PMC13276308",
"ISSN": "2589-0042",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.isci.2026.116320",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
11
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-72164-7 [code]
Multivariate genetic analysis reveals three distinct pathological dimensions in musculoskeletal disorders.
Journal: Nature communications
In common: data.table, tidyverse, genetics / omics, 2 references
[2] doi:10.1038/s41562-026-02486-5 [code]
Genome-wide association studies of infant and toddler temperament in European and multi-ancestry populations.
Journal: Nature human behaviour
In common: data.table, tidyverse, developmental, genetics / omics, 1 reference
[3] doi:10.1093/molbev/msag149 [code]
Multiple modes of selection underlie repeated and human-mediated adaptation in a formerly migratory fish.
Journal: Molecular biology and evolution
In common: data.table, tidyverse, other, 1 reference
[4] doi:10.1038/s42003-026-10131-0 [code]
Shared genetic architecture between the topology of brain white matter structural connectome and fluid intelligence.
Journal: Communications biology
In common: data.table, genetics / omics, 2 references
[5] doi:10.1038/s41514-026-00397-3 [code]
Nasal administration of Protollin enhances monocyte phagocytosis and decreases CD8&lt;sup&gt;+&lt;/sup&gt; T cell cytotoxicity in subjects with early Alzheimer's disease: a Phase 1 clinical trial.
Journal: npj aging
In common: data.table, tidyverse, 1 reference
[6] doi:10.3390/biomedicines14081677 [code]
Exploratory Genome and Transcriptome-Wide Association Analyses of Addiction-Related Phenotypes in a Twin Cohort.
Journal: Biomedicines
In common: data.table, tidyverse, genetics / omics, 1 reference
[7] doi:10.1038/s41588-026-02646-3 [code]
Co-expression-based models improve eQTL predictions for transcriptome-wide association studies and highlight new schizophrenia-associated genes.
Journal: Nature genetics
In common: data.table, tidyverse, genetics / omics, 1 reference
[8] doi:10.1038/s41467-026-73902-7 [code]
GWAS on short tandem repeats identifies genetic mechanisms in Alzheimer's disease.
Journal: Nature communications
In common: data.table, tidyverse, genetics / omics, 1 reference
[9] doi:10.1038/s41467-026-73428-y [code]
Regional heterogeneity in phenotypic and genetic associations between bone and brain in humans.
Journal: Nature communications
In common: data.table, tidyverse, genetics / omics, 1 reference
[10] doi:10.1186/s12967-026-08266-z [code]
Single-cell multi-omic integration analysis prioritizes druggable genes and reveals cell-type-specific causal effects in glioblastomagenesis.
Journal: Journal of translational medicine
In common: data.table, tidyverse, genetics / omics, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.