diff --git a/code/SoS/pecotmr_integration/twas_ctwas.ipynb b/code/SoS/pecotmr_integration/twas_ctwas.ipynb index 5375550fa..9bb1dc0b2 100644 --- a/code/SoS/pecotmr_integration/twas_ctwas.ipynb +++ b/code/SoS/pecotmr_integration/twas_ctwas.ipynb @@ -21,29 +21,13 @@ "source": [ "## Overview\n", "\n", - "A TWAS scan tests each gene's genetically predicted expression against a trait, but a\n", - "significant gene is not necessarily a causal one: nearby variants with direct effects on\n", - "the trait, and the predicted expression of neighbouring genes, are correlated with the\n", - "gene's own eQTLs and act as confounders. cTWAS addresses this by fine-mapping genes and\n", - "variants jointly within a region, so a gene is credited only for signal that its expression\n", - "explains beyond the surrounding variants and genes, and reports a posterior inclusion\n", - "probability rather than a p-value ([Zhao et al., 2024](https://doi.org/10.1038/s41588-023-01648-9)).\n", - "`finemapCtwasRegions`) and offers four workflows:\n", + "A TWAS scan tests each gene's genetically predicted expression, or another molecular trait, against a trait, but a significant gene is not necessarily a causal one: nearby variants with direct effects on the trait, and the predicted expression of neighbouring genes, are correlated with the gene's own xQTLs and act as confounders. cTWAS addresses this by fine-mapping genes and variants jointly within a region, so a gene is credited only for signal that its expression explains beyond the surrounding variants and genes, and reports a posterior inclusion probability rather than a p-value ([Zhao et al., 2024](https://doi.org/10.1038/s41588-023-01648-9)). This module follows the multi-group cTWAS formulation of Qian et al. (2024+), in which each molecular context forms its own group. Quantile TWAS extends the scan by testing genetic effects at different quantiles of the trait distribution. Mendelian randomization uses fine-mapped xQTL variants as instrumental variables (Zhang et al., 2020) to estimate a causal effect for candidate genes: it runs inside the same pecotmr pipeline when an xQTL fine-mapping result is supplied, aggregates multiple instruments by fixed-effect meta-analysis, and excludes results with severe exclusion-restriction violations. The unit of analysis is a single gene-trait pair.\n", "\n", - "This module provides software implementations for transcriptome-wide association analysis (TWAS), Quantile TWAS, and variant selection that yields sparse signals for cTWAS (causal TWAS) analysis following the multi-group cTWAS method of Qian et al. (2024+). It additionally performs Mendelian Randomization using fine-mapping instrumental variables (IV) as described in Zhang et al. (2020) for \"causal\" effect estimation and model validation. The unit of analysis is a single gene-trait pair.\n", + "TWAS and cTWAS read prediction weights stored as pecotmr `TwasWeights` objects. Weights from the legacy FunGen-xQTL release on Synapse use an older list layout and are converted first by the optional `twas_weights_conversion` workflow; weights already stored as `TwasWeights` skip this step. The conversion copies the weights and cross-validation metrics unchanged and writes the xQTL meta and data-type tables that the other workflows read. It is optional: skip it when the weights are already `TwasWeights`.\n", "\n", - "This notebook runs TWAS followed by cTWAS fine-mapping on a toy chr22 dataset: starting from pre-computed SuSiE-TWAS prediction weights for one gene, it tests each molecular context for association with a GWAS trait, keeps the imputable models selected by cross-validation, and then jointly fine-maps genes and SNPs within an LD block to identify which signals are likely to be directly causal rather than driven by correlation. Quantile TWAS extends traditional TWAS by testing genetic effects at different quantiles of the trait distribution.\n", + "The module offers four workflows, in the order they run: `twas_weights_conversion` (optional) converts legacy weights, `twas` tests each gene and molecular context for association using the models that pass the cross-validation cutoffs, `ctwas` fine-maps genes and variants jointly in three chained steps, and `quantile_twas` is an alternative scan across trait quantiles.\n", "\n", - "This procedure is a continuation of the SuSiE-TWAS workflow: it assumes that xQTL fine-mapping has been performed and molecular-trait prediction weights pre-computed (to be used for TWAS). Cross-validation of TWAS weights is optional but highly recommended.\n", - "\n", - "**Prerequisites:** TWAS weights (`protocol_example.twas_weights.rds`) from `mnm_regression`, plus GWAS summary statistics and an LD matrix for the region of interest.\n", - "\n", - "**Timing** (toy chr22 dataset): `twas` ~50 sec; `ctwas` ~42 sec.\n", - "\n", - "MR can be run for candidate genes afterwards, limited to genes with cTWAS significance and a strong instrumental variable, using fine-mapped xQTL with GWAS data; multiple IVs aggregate by fixed-effect meta-analysis, and results with severe exclusion-restriction violations are excluded.\n", - "\n", - "**When to run it.** After you have TWAS weights (from `mnm_regression`) and GWAS summary\n", - "statistics. cTWAS consumes weights, it does not produce them." + "**When to run it.** After xQTL fine-mapping and TWAS weight training (`mnm_regression`), with GWAS summary statistics and an LD reference in hand. Weights from the legacy FunGen-xQTL release need `twas_weights_conversion` first; weights already in pecotmr `TwasWeights` format do not. Cross-validated weights are strongly recommended. cTWAS consumes weights; it does not produce them." ] }, { @@ -63,12 +47,12 @@ "protocol_example_twas_chr22\t22\tprotocol_example.twas.gwas_sumstats.chr22.tsv.gz\t200000\n", "```\n", "\n", - "- `--xqtl_meta_data` **`tests/fixtures/twas/protocol_example.twas.xqtl_meta.tsv`**\n", - "(one row per gene or region with its weight file. Columns `#chr`, `start`, `end`, `region_id`, `TSS`, `original_data`, `contexts`.)\n", + "- `--xqtl_meta_data` **`output/twas_weights_conversion/protocol_example.protocol_example.xqtl_meta.tsv`**\n", + "(one row per gene with its `TwasWeights` file, as written by `twas_weights_conversion`; for weights already in that format, a table of the same layout pointing at them. Columns `#chr`, `start`, `end`, `region_id`, `TSS`, `original_data`, `contexts`.)\n", "\n", "```\n", "#chr\tstart\tend\tregion_id\tTSS\toriginal_data\tcontexts\n", - "chr22\t10000000\t19000000\tENSG00000130538\t15528191\ttests/fixtures/twas/protocol_example.twas.reshaped_toy.chr22_ENSG00000130538.univariate_twas_weights.rds\tbulk_rnaseq\n", + "chr22\t0\t18960000\tENSG00000130538\t10685239\t/output/twas_weights_conversion/protocol_example/protocol_example.ENSG00000130538.twas_weights.s4.rds\tbulk_rnaseq\n", "```\n", "\n", "- `--ld_meta_data` **`tests/fixtures/ld_reference/ld_meta_file.tsv`**\n", @@ -88,14 +72,22 @@ "chr22\t10000000\t19000000\n", "```\n", "\n", - "- `--xqtl_type_table` **`tests/fixtures/twas/protocol_example.twas.data_type_table.txt`**\n", - "(maps each xQTL context to its modality, used when grouping contexts for cTWAS)\n", + "- `--xqtl_type_table` **`output/twas_weights_conversion/protocol_example.xqtl_type_table.tsv`**\n", + "(maps each xQTL context to its modality, used when grouping contexts for cTWAS; written by `twas_weights_conversion`, or supplied directly)\n", "\n", "```\n", "context\ttype\n", "bulk_rnaseq\teQTL\n", "```\n", "\n", + "- `--legacy_weights` **`tests/fixtures/twas/protocol_example.twas.legacy_weights.tsv`** (`twas_weights_conversion` only)\n", + "(legacy FunGen-xQTL weight files to convert, one per row. Columns `study`, `data_type`, `path`; rows that share a study and a gene are merged into one weight file.)\n", + "\n", + "```\n", + "study\tdata_type\tpath\n", + "protocol_example\teQTL\ttests/fixtures/twas/protocol_example.twas.reshaped_toy.chr22_ENSG00000130538.univariate_twas_weights.rds\n", + "```\n", + "\n", "- `--cwd output/twas` (working directory all outputs are written under)\n", "- `--name protocol_example` (tag used in every output filename)\n", "- `--region-name chr22_10000000_19000000` (optional; restrict to named regions)\n", @@ -104,6 +96,7 @@ "\n", "- `--ld_reference_sample_size` (sample size of the LD reference panel, used to scale the LD matrices)\n", "- `--prior_var_structure` (how the cTWAS prior variance is shared across groups)\n", + "- `--gwas_study ` (`ctwas` only: the one GWAS study to fine-map, required when `--gwas_meta_data` lists several; run each study with its own `--cwd`)\n", "- `--rsq_cutoff` (minimum cross-validation r-squared for a weight model to be used)\n", "- `--rsq_pval_cutoff` (maximum cross-validation p-value for a weight model to be used)\n", "- `--twas_weight_cutoff` (drop weights below this magnitude before assembling cTWAS inputs)\n", @@ -135,6 +128,24 @@ "source": [ "## Output\n", "\n", + "- **`output/twas_weights_conversion//..twas_weights.s4.rds`** (`twas_weights_conversion`)\n", + "(one pecotmr `TwasWeights` object per study and gene; one row per context and method, weights copied unchanged from the legacy file)\n", + "- **`output/twas_weights_conversion/..xqtl_meta.tsv`** (`twas_weights_conversion`)\n", + "(`--xqtl_meta_data` input of `twas` and `ctwas`, one table per study)\n", + "\n", + "```\n", + "#chr\tstart\tend\tregion_id\tTSS\toriginal_data\tcontexts\n", + "chr22\t0\t18960000\tENSG00000130538\t10685239\t/output/twas_weights_conversion/protocol_example/protocol_example.ENSG00000130538.twas_weights.s4.rds\tbulk_rnaseq\n", + "```\n", + "\n", + "- **`output/twas_weights_conversion/.xqtl_type_table.tsv`** (`twas_weights_conversion`)\n", + "(`--xqtl_type_table` input: the data type of every converted context)\n", + "\n", + "```\n", + "context\ttype\n", + "bulk_rnaseq\teQTL\n", + "```\n", + "\n", "- **`output/twas/twas/...twas.rds`** (`twas`)\n", "(per-gene TWAS result, one file per gene in the region)\n", "- **`output/twas/twas/..twas.tsv.gz`** (`twas`, side file)\n", @@ -173,7 +184,46 @@ "source": [ "## Minimal Working Example\n", "\n", - "The three methods below are independent: run whichever you need. All use the chr22 example data in this repository." + "The three methods below are independent: run whichever you need. All use the chr22 example data in this repository.\n", + "\n", + "`twas_weights_conversion`, below, is only needed for weights in the legacy FunGen-xQTL format, as the example weights here are; its tables are the `--xqtl_meta_data` and `--xqtl_type_table` inputs of the three methods. Weights already stored as pecotmr `TwasWeights` skip it and go straight to the three methods with a table of the same layout pointing at them." + ] + }, + { + "cell_type": "markdown", + "id": "1f7f53b5", + "metadata": { + "kernel": "SoS" + }, + "source": [ + "### Convert legacy TWAS weights\n", + "\n", + "`twas_weights_conversion` converts legacy FunGen-xQTL weight files into pecotmr `TwasWeights`, one per study and gene, and writes the xQTL meta and data-type tables that the other workflows read." + ] + }, + { + "cell_type": "markdown", + "id": "8dce9984", + "metadata": { + "kernel": "SoS" + }, + "source": [ + "**Timing**: ~20 sec (on toy dataset)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f6b00db", + "metadata": { + "kernel": "SoS" + }, + "outputs": [], + "source": [ + "sos run pipeline/twas_ctwas.ipynb twas_weights_conversion \\\n", + " --cwd output \\\n", + " --name protocol_example \\\n", + " --legacy_weights tests/fixtures/twas/protocol_example.twas.legacy_weights.tsv" ] }, { @@ -195,7 +245,7 @@ "kernel": "SoS" }, "source": [ - "**Timing**: TBD (on toy dataset)" + "**Timing**: ~25 sec (on toy dataset)" ] }, { @@ -210,11 +260,11 @@ "sos run pipeline/twas_ctwas.ipynb twas \\\n", " --cwd output --name protocol_example \\\n", " --gwas_meta_data tests/fixtures/twas/protocol_example.twas.gwas_meta.tsv \\\n", - " --xqtl_meta_data tests/fixtures/twas/protocol_example.twas.xqtl_meta.tsv \\\n", + " --xqtl_meta_data output/twas_weights_conversion/protocol_example.protocol_example.xqtl_meta.tsv \\\n", " --ld_meta_data tests/fixtures/ld_reference/ld_meta_file.tsv \\\n", " --ld_reference_sample_size 17000 \\\n", " --regions tests/fixtures/twas/protocol_example.twas.LD_blocks.chr22.bed \\\n", - " --xqtl_type_table tests/fixtures/twas/protocol_example.twas.data_type_table.txt \\\n", + " --xqtl_type_table output/twas_weights_conversion/protocol_example.xqtl_type_table.tsv \\\n", " --rsq_pval_cutoff 0.05 --rsq_cutoff 0.01 \\\n", " --region-name chr22_10000000_19000000" ] @@ -268,7 +318,7 @@ " --cwd output --name protocol_example \\\n", " --thin 1 --prior_var_structure shared_all \\\n", " --gwas_meta_data tests/fixtures/twas/protocol_example.twas.gwas_meta.tsv \\\n", - " --xqtl_meta_data tests/fixtures/twas/protocol_example.twas.xqtl_meta.tsv \\\n", + " --xqtl_meta_data output/twas_weights_conversion/protocol_example.protocol_example.xqtl_meta.tsv \\\n", " --ld_meta_data tests/fixtures/ld_reference/ld_meta_file.tsv \\\n", " --regions tests/fixtures/twas/protocol_example.twas.LD_blocks.chr22.bed \\\n", " --twas_weight_cutoff 0 \\\n", @@ -309,7 +359,7 @@ " --prior_var_structure shared_all \\\n", " --cwd output --name protocol_example \\\n", " --gwas_meta_data tests/fixtures/twas/protocol_example.twas.gwas_meta.tsv \\\n", - " --xqtl_meta_data tests/fixtures/twas/protocol_example.twas.xqtl_meta.tsv \\\n", + " --xqtl_meta_data output/twas_weights_conversion/protocol_example.protocol_example.xqtl_meta.tsv \\\n", " --ld_meta_data tests/fixtures/ld_reference/ld_meta_file.tsv \\\n", " --regions tests/fixtures/twas/protocol_example.twas.LD_blocks.chr22.bed" ] @@ -348,7 +398,7 @@ " --prior_var_structure shared_all \\\n", " --cwd output --name protocol_example \\\n", " --gwas_meta_data tests/fixtures/twas/protocol_example.twas.gwas_meta.tsv \\\n", - " --xqtl_meta_data tests/fixtures/twas/protocol_example.twas.xqtl_meta.tsv \\\n", + " --xqtl_meta_data output/twas_weights_conversion/protocol_example.protocol_example.xqtl_meta.tsv \\\n", " --ld_meta_data tests/fixtures/ld_reference/ld_meta_file.tsv \\\n", " --regions tests/fixtures/twas/protocol_example.twas.LD_blocks.chr22.bed \\\n", " --region-name chr22_10000000_19000000" @@ -387,11 +437,11 @@ "sos run pipeline/twas_ctwas.ipynb quantile_twas \\\n", " --cwd output --name protocol_example \\\n", " --gwas_meta_data tests/fixtures/twas/protocol_example.twas.gwas_meta.tsv \\\n", - " --xqtl_meta_data tests/fixtures/twas/protocol_example.twas.xqtl_meta.tsv \\\n", + " --xqtl_meta_data output/twas_weights_conversion/protocol_example.protocol_example.xqtl_meta.tsv \\\n", " --ld_meta_data tests/fixtures/ld_reference/ld_meta_file.tsv \\\n", " --ld_reference_sample_size 17000 \\\n", " --regions tests/fixtures/twas/protocol_example.twas.LD_blocks.chr22.bed \\\n", - " --xqtl_type_table tests/fixtures/twas/protocol_example.twas.data_type_table.txt \\\n", + " --xqtl_type_table output/twas_weights_conversion/protocol_example.xqtl_type_table.tsv \\\n", " --region-name chr22_10000000_19000000" ] }, @@ -433,6 +483,7 @@ " workflow_options: Double-hyphen workflow-specific parameters\n", "\n", "Workflows:\n", + " twas_weights_conversion\n", " get_analysis_regions\n", " twas\n", " ctwas\n", @@ -469,6 +520,22 @@ " specific analysis\n", "\n", "Sections\n", + " twas_weights_conversion_1:\n", + " Workflow Options:\n", + " --legacy-weights . (as path)\n", + " Convert legacy FunGen-xQTL TWAS weight RDS files (the\n", + " list layout released on Synapse: gene ->\n", + " \"_\" -> twas_weights / twas_cv_result /\n", + " region_info) into pecotmr TwasWeights, one RDS per\n", + " (study, gene), with twas_weights_conversion.R. Weights\n", + " and CV metrics are copied unchanged; no LD sketch is\n", + " attached. --legacy-weights is a TSV with columns study,\n", + " data_type, path (one legacy RDS per row). Rows sharing a\n", + " study and a gene (e.g. the univariate and multicontext\n", + " files of one gene) are merged into one TwasWeights,\n", + " because the twas step reads one weight file per gene.\n", + " The gene is the ENSG ID in the legacy file name.\n", + " twas_weights_conversion_2:\n", " get_analysis_regions:\n", " twas:\n", " Workflow Options:\n", @@ -621,6 +688,87 @@ "\n" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "12c79b7f", + "metadata": { + "kernel": "SoS" + }, + "outputs": [], + "source": [ + "[twas_weights_conversion_1]\n", + "# Convert legacy FunGen-xQTL TWAS weight RDS files (the list layout released on\n", + "# Synapse: gene -> \"_\" -> twas_weights / twas_cv_result /\n", + "# region_info) into pecotmr TwasWeights, one RDS per (study, gene), with\n", + "# twas_weights_conversion.R. Weights and CV metrics are copied unchanged; no LD\n", + "# sketch is attached. --legacy-weights is a TSV with columns study, data_type,\n", + "# path (one legacy RDS per row). Rows sharing a study and a gene (e.g. the\n", + "# univariate and multicontext files of one gene) are merged into one\n", + "# TwasWeights, because the twas step reads one weight file per gene. The gene\n", + "# is the ENSG ID in the legacy file name.\n", + "parameter: legacy_weights = path()\n", + "import csv, re\n", + "_legacy_rows = list(csv.DictReader(open(legacy_weights), delimiter = '\\t'))\n", + "stop_if(len(_legacy_rows) == 0, f\"No rows in {legacy_weights}.\")\n", + "_conv_units = {}\n", + "for _r in _legacy_rows:\n", + " _g = re.search(r'ENSG[0-9]+', path(_r['path']).name)\n", + " stop_if(_g is None, f\"No ENSG gene ID in legacy file name: {_r['path']}\")\n", + " _conv_units.setdefault((_r['study'], _g.group(0)), []).append(_r['path'])\n", + "_conv_keys = sorted(_conv_units)\n", + "_conv_types = {_r['study']: _r['data_type'] for _r in _legacy_rows}\n", + "input: [f for k in _conv_keys for f in _conv_units[k]], group_by = lambda x: [_conv_units[k] for k in _conv_keys]\n", + "output: f'{cwd:a}/twas_weights_conversion/{_conv_keys[_index][0]}/{_conv_keys[_index][0]}.{_conv_keys[_index][1]}.twas_weights.s4.rds'\n", + "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\n", + "bash: expand = '${ }', stdout = f\"{_output:n}.stdout\", stderr = f\"{_output:n}.stderr\", container = container\n", + " Rscript ${modular_script_dir}/pecotmr_integration/twas_weights_conversion.R \\\n", + " --legacy \"${\",\".join([str(x) for x in _input])}\" \\\n", + " --study \"${_conv_keys[_index][0]}\" \\\n", + " --data-type \"${_conv_types[_conv_keys[_index][0]]}\" \\\n", + " --output \"${_output}\" \\\n", + " --meta-row \"${_output:n}.meta_row.tsv\" \\\n", + " --type-rows \"${_output:n}.type_rows.tsv\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "094ba82e", + "metadata": { + "kernel": "SoS" + }, + "outputs": [], + "source": [ + "[twas_weights_conversion_2]\n", + "# Collect the per-gene rows written by step 1 into one xqtl_meta table per study\n", + "# ({name}..xqtl_meta.tsv, the --xqtl-meta-data input of twas and ctwas)\n", + "# and one context -> data type table ({name}.xqtl_type_table.tsv, the\n", + "# --xqtl-type-table input). The meta tables are per study because twas uses\n", + "# the first weight file it finds for each gene: run twas once per study.\n", + "input: group_by = 'all'\n", + "_conv_studies = sorted(set(path(x).parent.name for x in _input))\n", + "output: [f'{cwd:a}/twas_weights_conversion/{name}.{s}.xqtl_meta.tsv' for s in _conv_studies] + [f'{cwd:a}/twas_weights_conversion/{name}.xqtl_type_table.tsv']\n", + "_ctx_type = {}\n", + "for _s, _out in zip(_conv_studies, _output[:-1]):\n", + " _rows = []\n", + " for _x in _input:\n", + " if path(_x).parent.name != _s:\n", + " continue\n", + " _rows += [l.rstrip('\\n').split('\\t') for l in open(f'{_x:n}.meta_row.tsv') if l.strip()]\n", + " for l in open(f'{_x:n}.type_rows.tsv'):\n", + " if l.strip():\n", + " _c, _t = l.rstrip('\\n').split('\\t')\n", + " _ctx_type[_c] = _t\n", + " _rows.sort(key = lambda r: (r[0], int(r[1]), r[3]))\n", + " with open(_out, 'w') as _f:\n", + " _f.write('#chr\\tstart\\tend\\tregion_id\\tTSS\\toriginal_data\\tcontexts\\n')\n", + " _f.writelines('\\t'.join(r) + '\\n' for r in _rows)\n", + "with open(_output[-1], 'w') as _f:\n", + " _f.write('context\\ttype\\n')\n", + " _f.writelines(f'{c}\\t{t}\\n' for c, t in sorted(_ctx_type.items()))" + ] + }, { "cell_type": "code", "execution_count": null, @@ -789,6 +937,10 @@ "parameter: max_num_variants = \"Inf\"\n", "parameter: cs_min_cor = 0.0\n", "parameter: min_pip_cutoff = 0.0\n", + "# TWAS method fed to cTWAS (assembleCtwasInputs method=). Empty: pecotmr picks\n", + "# \"ensemble\" if present, else the sole method; weights carrying several methods and\n", + "# no ensemble (e.g. converted legacy weights) need one, e.g. --ctwas-method susie.\n", + "parameter: ctwas_method = \"\"\n", "# declared for CLI stability; consumed downstream / not by assembleCtwasInputs\n", "parameter: thin = 1.0\n", "parameter: maxSNP = 20000\n", @@ -809,7 +961,24 @@ "_chrom_jobs = [j for j in jobs if j[\"chrom\"] == ctwas_chrom]\n", "ctwas_studies = _chrom_jobs[0][\"gwas_studies\"].split(',') if _chrom_jobs and _chrom_jobs[0][\"gwas_studies\"] else []\n", "ctwas_gwas_files = _chrom_jobs[0][\"gwas_files\"] if _chrom_jobs else \"\"\n", + "# assembleCtwasInputs() takes one GWAS per LD block, so cTWAS runs one study at\n", + "# a time: --gwas-study picks it when --gwas-meta-data lists several. The\n", + "# per-block GWAS files are named by region only, so give each study its own --cwd.\n", + "if gwas_study:\n", + " _keep = [i for i, s in enumerate(ctwas_studies) if s in gwas_study]\n", + " ctwas_gwas_files = \",\".join([ctwas_gwas_files.split(\",\")[i] for i in _keep])\n", + " ctwas_studies = [ctwas_studies[i] for i in _keep]\n", "stop_if(len(ctwas_studies) == 0, f\"No GWAS study covers {ctwas_chrom}.\")\n", + "fail_if(len(ctwas_studies) > 1, f\"cTWAS takes one GWAS study at a time; {ctwas_chrom} has {ctwas_studies}. Pass --gwas-study.\")\n", + "# ctwasPipeline needs the LD panel as a real genotype file: a block built from the\n", + "# chromosome->path mapping of --ld-meta-data keeps only a \"\" placeholder.\n", + "# Resolve this chromosome's genotype prefix and pass it as --ld-sketch when unique.\n", + "import csv as _csv, os as _os\n", + "_ld_rows = list(_csv.DictReader(open(ld_meta_data), delimiter = \"\\t\"))\n", + "_ld_chr = [x for x in (\"#chr\", \"#chrom\", \"chr\", \"chrom\") if _ld_rows and x in _ld_rows[0]]\n", + "_ld_dir = _os.path.dirname(_os.path.abspath(ld_meta_data))\n", + "_ld_paths = sorted(set(_os.path.join(_ld_dir, r[\"path\"].split(\",\")[0]) for r in _ld_rows if _ld_chr and r[_ld_chr[0]].replace(\"chr\", \"\") == ctwas_chrom.replace(\"chr\", \"\")))\n", + "ctwas_ld_sketch = _ld_paths[0] if len(_ld_paths) == 1 and any(_os.path.exists(_ld_paths[0] + e) for e in (\".pgen\", \".bed\")) else \"\"\n", "output: f\"{cwd:a}/ctwas/{name}.ctwas_inputs.rds\"\n", "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f\"{step_name}_{_output:bn}\"\n", "bash: expand = '${ }', stdout = f\"{_output:n}.stdout\", stderr = f\"{_output:n}.stderr\", container = container\n", @@ -834,13 +1003,14 @@ " --gwas-tsv \"$gwas_tsvs\" \\\n", " --ld-block \"$region\" \\\n", " --ld-meta \"${ld_meta_data}\" \\\n", - " --output \"$gwas_rds\"\n", + " ${f'--ld-sketch \"{ctwas_ld_sketch}\" ' if ctwas_ld_sketch else \"\"}--output \"$gwas_rds\"\n", " done\n", "\n", " # (3) assemble cTWAS inputs (pecotmr places each gene into its home block).\n", " Rscript ${modular_script_dir}/pecotmr_integration/ctwas_assemble.R \\\n", " --manifest \"$manifest\" \\\n", " --twas-weights \"${\",\".join([str(w) for w in ctwas_weights])}\" \\\n", + " ${f\"--method {ctwas_method} \" if ctwas_method else \"\"}\\\n", " --twas-weight-cutoff ${twas_weight_cutoff} \\\n", " --cs-min-cor ${cs_min_cor} \\\n", " --min-pip-cutoff ${min_pip_cutoff} \\\n", diff --git a/code/script/pecotmr_integration/ctwas_assemble.R b/code/script/pecotmr_integration/ctwas_assemble.R index 709408e39..427f0a184 100644 --- a/code/script/pecotmr_integration/ctwas_assemble.R +++ b/code/script/pecotmr_integration/ctwas_assemble.R @@ -97,12 +97,33 @@ if (nrow(manifest) < 2L) # region_id (a per-block construct call keys on the seqname, which is the same # for every block on a chromosome) and combined. combineGwasSumStats() unions # the per-block LD panels and concatenates the per-element QC audit. -gwasParts <- vector("list", nrow(manifest)) +# +# A study with no variants in a block (QC gaps, centromeres, chromosome ends) +# leaves an empty element: it carries no genome build and an empty range, so +# combineGwasSumStats() rejects it ("must share one genome build (got NA)" or +# "(study, range) must be unique"). Drop empty elements, and blocks left with +# none. Subsetting drops the ldSketch, qcInfo and metadata slots; restore them. +gwasParts <- list() for (i in seq_len(nrow(manifest))) { part <- readRDS(manifest$gwas_sumstats_rds[[i]]) + keep <- lengths(part) > 0L + if (!any(keep)) { + message("Skipping LD block ", manifest$region_id[[i]], ": no GWAS variants.") + next + } + if (!all(keep)) { + message("LD block ", manifest$region_id[[i]], ": dropping ", sum(!keep), + " study element(s) with no variants.") + kept <- part[keep] + for (s in c("ldSketch", "qcInfo", "metadata")) + methods::slot(kept, s) <- methods::slot(part, s) + part <- kept + } S4Vectors::mcols(part)$blockId <- manifest$region_id[[i]] - gwasParts[[i]] <- part + gwasParts[[length(gwasParts) + 1L]] <- part } +if (length(gwasParts) < 2L) + stop("Fewer than two LD blocks have GWAS variants; cTWAS's EM needs multi-block context.") gwasSumStats <- combineGwasSumStats(gwasParts) # FLAT weight source: an unnamed list of per-gene weight objects. assembleCtwasInputs diff --git a/code/script/pecotmr_integration/twas_weights_conversion.R b/code/script/pecotmr_integration/twas_weights_conversion.R new file mode 100644 index 000000000..ce008d2a7 --- /dev/null +++ b/code/script/pecotmr_integration/twas_weights_conversion.R @@ -0,0 +1,91 @@ +#!/usr/bin/env Rscript +# twas_weights_conversion.R +# +# Convert legacy FunGen-xQTL TWAS weight RDS files (the frozen Synapse release; +# list layout gene -> "_" -> twas_weights / twas_cv_result / +# region_info) into ONE pecotmr (>= 0.8.2) TwasWeights object per gene, and +# write the matching xqtl_meta row +# (#chr start end region_id TSS original_data contexts). +# +# Several legacy files for the same gene (e.g. univariate + multicontext) are +# merged into one TwasWeights, because twas_ctwas.ipynb uses one weight file +# per gene. No LD sketch is attached (ldSketch = NULL), as in the protocol +# fixture protocol_example.twas.reshaped_toy.*.twas_weights.s4.rds. +# +# Usage: +# Rscript twas_weights_conversion.R \ +# --legacy a.univariate_twas_weights.rds,a.multicontext_twas_weights.rds \ +# --study MSBB_eQTL --output gene.twas_weights.s4.rds --meta-row gene.meta.tsv +# [--data-type eQTL --type-rows gene.types.tsv] + +suppressPackageStartupMessages({ library(argparser); library(pecotmr) }) +p <- arg_parser("Legacy TWAS weights -> pecotmr TwasWeights") +p <- add_argument(p, "--legacy", help = "Comma-separated legacy RDS files for ONE gene") +p <- add_argument(p, "--study", help = "Study label, e.g. MSBB_eQTL") +p <- add_argument(p, "--output", help = "Output TwasWeights RDS") +p <- add_argument(p, "--meta-row", help = "Output one-row xqtl_meta TSV (no header)") +p <- add_argument(p, "--data-type", default = "", help = "Optional data type of this study (eQTL, pQTL, sQTL, ...)") +p <- add_argument(p, "--type-rows", default = "", help = "Optional output: contextdata_type rows (no header); needs --data-type") +a <- parse_args(p) + +files <- trimws(strsplit(a$legacy, ",")[[1]]) +S <- C <- Tr <- M <- character(0); E <- list(); ctxs <- character(0); info <- NULL +# traitPos: the gene TSS for every row; cTWAS places each gene into its LD block by it +TPc <- character(0); TPp <- integer(0) + +# legacy CV performance (1-row matrix: corr rsq adj_rsq pval RMSE MAE) +# -> list(metrics = named numeric), the layout used by the protocol fixture +cvMetrics <- function(perf) { + if (is.null(perf)) return(NULL) + v <- unlist(as.data.frame(perf)[1, , drop = TRUE]) + list(metrics = setNames(as.numeric(v), names(v))) +} + +for (f in files) { + lg <- readRDS(f) + for (gene in names(lg)) for (cn in names(lg[[gene]])) { + co <- lg[[gene]][[cn]] + if (!is.list(co) || is.null(co$twas_weights)) next + if (is.null(info)) info <- co$region_info + ctx <- sub(paste0("[_:]", gene, "$"), "", cn) # sQTL names end in ":" + perf <- co$twas_cv_result$performance + names(perf) <- sub("_performance$", "", names(perf)) + for (wn in names(co$twas_weights)) { + w <- co$twas_weights[[wn]] + vids <- if (is.matrix(w)) rownames(w) else names(w) + if (is.null(vids)) vids <- co$variant_names + wv <- if (is.matrix(w)) as.numeric(w[, 1]) else as.numeric(w) + if (all(is.na(wv)) || all(wv == 0, na.rm = TRUE)) next + m <- sub("_weights$", "", wn) + fits <- if (m %in% c("susie", "susie_inf", "susie_ash")) co$susie_weights_intermediate else NULL + E[[length(E) + 1L]] <- twasWeightsRow(variantIds = vids, weights = wv, + fits = fits, cvResult = cvMetrics(perf[[m]])) + S <- c(S, a$study); C <- c(C, ctx); Tr <- c(Tr, gene); M <- c(M, m) + TPc <- c(TPc, paste0("chr", sub("^chr", "", co$region_info$region_coord$chrom))) + TPp <- c(TPp, as.integer(co$region_info$region_coord$start)) + } + ctxs <- union(ctxs, ctx) + } +} +if (!length(E)) stop("No usable (non-zero) weights in: ", a$legacy) + +tw <- TwasWeights(study = S, context = C, trait = Tr, method = M, entry = E, + traitPos = GenomicRanges::GRanges(TPc, IRanges::IRanges(TPp, width = 1L))) +dir.create(dirname(a$output), showWarnings = FALSE, recursive = TRUE) +saveRDS(tw, a$output) + +g <- info$grange; rc <- info$region_coord +# Released weights store the TSS as a 1-bp region_coord; warn if it is a range +if (rc$end - rc$start > 1) + warning("region_coord spans ", rc$start, "-", rc$end, " (not a single TSS); using its start as TSS") +row <- data.frame(paste0("chr", sub("^chr", "", g$chrom)), g$start, g$end, + unique(Tr)[1], rc$start, normalizePath(a$output), + paste(ctxs, collapse = ",")) +write.table(row, a$meta_row, sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE) +if (nzchar(a$type_rows)) { + if (!nzchar(a$data_type)) stop("--type-rows needs --data-type") + write.table(data.frame(ctxs, a$data_type), a$type_rows, sep = "\t", quote = FALSE, + row.names = FALSE, col.names = FALSE) +} +cat(sprintf("Wrote %d weight rows (%d contexts) for %s -> %s\n", + length(E), length(ctxs), unique(Tr)[1], a$output)) diff --git a/tests/fixtures/twas/protocol_example.twas.legacy_weights.tsv b/tests/fixtures/twas/protocol_example.twas.legacy_weights.tsv new file mode 100644 index 000000000..78ff655f6 --- /dev/null +++ b/tests/fixtures/twas/protocol_example.twas.legacy_weights.tsv @@ -0,0 +1,2 @@ +study data_type path +protocol_example eQTL tests/fixtures/twas/protocol_example.twas.reshaped_toy.chr22_ENSG00000130538.univariate_twas_weights.rds