From 35e93fdd87e36befce15a18636d2813a369bcde2 Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 8 Jul 2026 03:32:26 -0400 Subject: [PATCH 01/12] fix(genotype-preprocessing): replace unused name placeholder with required cohort/context parameters GWAS_QC, PCA, genotype_formatting, genotype_alignment, and TensorQTL had a 'parameter: name = ""' (or similarly unused) placeholder that was never wired to the documented --name example, so the MWE commands silently used an empty run tag. Replaced with explicit required parameters (--cohort, and --cohort/--context/--modality for TensorQTL) that are actually used to build the output name, matching the protocol_example naming convention. --- .../TensorQTL/TensorQTL.ipynb | 49 ++++++++++--------- .../data_preprocessing/genotype/GWAS_QC.ipynb | 11 +++-- .../SoS/data_preprocessing/genotype/PCA.ipynb | 40 ++++++++------- .../genotype/genotype_formatting.ipynb | 6 ++- code/SoS/misc/genotype_alignment.ipynb | 7 +-- 5 files changed, 63 insertions(+), 50 deletions(-) diff --git a/code/SoS/association_scan/TensorQTL/TensorQTL.ipynb b/code/SoS/association_scan/TensorQTL/TensorQTL.ipynb index 8ba088f90..48b60818a 100644 --- a/code/SoS/association_scan/TensorQTL/TensorQTL.ipynb +++ b/code/SoS/association_scan/TensorQTL/TensorQTL.ipynb @@ -56,7 +56,7 @@ "data": { "text/html": [ "\n", - "\n", + "\n", "\n", "\t\n", "\t\n", @@ -67,7 +67,7 @@ "
A data.table: 1 × 2A data.table: 1 \u00d7 2
#id#dir
<int><chr>
\n" ], "text/latex": [ - "A data.table: 1 × 2\n", + "A data.table: 1 \u00d7 2\n", "\\begin{tabular}{ll}\n", " \\#id & \\#dir\\\\\n", " & \\\\\n", @@ -77,7 +77,7 @@ ], "text/markdown": [ "\n", - "A data.table: 1 × 2\n", + "A data.table: 1 \u00d7 2\n", "\n", "| #id <int> | #dir <chr> |\n", "|---|---|\n", @@ -96,7 +96,7 @@ "data": { "text/html": [ "\n", - "\n", + "\n", "\n", "\t\n", "\t\n", @@ -111,7 +111,7 @@ "
A data.table: 5 × 8A data.table: 5 \u00d7 8
#chrstartendIDSAMPLE_001SAMPLE_002SAMPLE_003SAMPLE_004
<chr><int><int><chr><dbl><dbl><dbl><dbl>
\n" ], "text/latex": [ - "A data.table: 5 × 8\n", + "A data.table: 5 \u00d7 8\n", "\\begin{tabular}{llllllll}\n", " \\#chr & start & end & ID & SAMPLE\\_001 & SAMPLE\\_002 & SAMPLE\\_003 & SAMPLE\\_004\\\\\n", " & & & & & & & \\\\\n", @@ -125,7 +125,7 @@ ], "text/markdown": [ "\n", - "A data.table: 5 × 8\n", + "A data.table: 5 \u00d7 8\n", "\n", "| #chr <chr> | start <int> | end <int> | ID <chr> | SAMPLE_001 <dbl> | SAMPLE_002 <dbl> | SAMPLE_003 <dbl> | SAMPLE_004 <dbl> |\n", "|---|---|---|---|---|---|---|---|\n", @@ -192,7 +192,7 @@ "data": { "text/html": [ "\n", - "\n", + "\n", "\n", "\t\n", "\t\n", @@ -208,7 +208,7 @@ "
A data.table: 6 × 2A data.table: 6 \u00d7 2
#id#path
<int><chr>
\n" ], "text/latex": [ - "A data.table: 6 × 2\n", + "A data.table: 6 \u00d7 2\n", "\\begin{tabular}{ll}\n", " \\#id & \\#path\\\\\n", " & \\\\\n", @@ -223,7 +223,7 @@ ], "text/markdown": [ "\n", - "A data.table: 6 × 2\n", + "A data.table: 6 \u00d7 2\n", "\n", "| #id <int> | #path <chr> |\n", "|---|---|\n", @@ -271,7 +271,7 @@ "data": { "text/html": [ "\n", - "\n", + "\n", "\n", "\t\n", "\n", @@ -285,7 +285,7 @@ "
A matrix: 5 × 5 of type intA matrix: 5 \u00d7 5 of type int
SAMPLE_001SAMPLE_002SAMPLE_003SAMPLE_004SAMPLE_005
\n" ], "text/latex": [ - "A matrix: 5 × 5 of type int\n", + "A matrix: 5 \u00d7 5 of type int\n", "\\begin{tabular}{r|lllll}\n", " & SAMPLE\\_001 & SAMPLE\\_002 & SAMPLE\\_003 & SAMPLE\\_004 & SAMPLE\\_005\\\\\n", "\\hline\n", @@ -298,7 +298,7 @@ ], "text/markdown": [ "\n", - "A matrix: 5 × 5 of type int\n", + "A matrix: 5 \u00d7 5 of type int\n", "\n", "| | SAMPLE_001 | SAMPLE_002 | SAMPLE_003 | SAMPLE_004 | SAMPLE_005 |\n", "|---|---|---|---|---|---|\n", @@ -356,7 +356,7 @@ "data": { "text/html": [ "\n", - "\n", + "\n", "\n", "\t\n", "\t\n", @@ -394,7 +394,7 @@ "
A data.table: 28 × 5A data.table: 28 \u00d7 5
#idSAMPLE_001SAMPLE_002SAMPLE_003SAMPLE_004
<chr><dbl><dbl><dbl><dbl>
\n" ], "text/latex": [ - "A data.table: 28 × 5\n", + "A data.table: 28 \u00d7 5\n", "\\begin{tabular}{lllll}\n", " \\#id & SAMPLE\\_001 & SAMPLE\\_002 & SAMPLE\\_003 & SAMPLE\\_004\\\\\n", " & & & & \\\\\n", @@ -431,7 +431,7 @@ ], "text/markdown": [ "\n", - "A data.table: 28 × 5\n", + "A data.table: 28 \u00d7 5\n", "\n", "| #id <chr> | SAMPLE_001 <dbl> | SAMPLE_002 <dbl> | SAMPLE_003 <dbl> | SAMPLE_004 <dbl> |\n", "|---|---|---|---|---|\n", @@ -559,7 +559,7 @@ "data": { "text/html": [ "\n", - "\n", + "\n", "\n", "\t\n", "\t\n", @@ -575,7 +575,7 @@ "
A data.table: 6 × 4A data.table: 6 \u00d7 4
#chrstartendgene_id
<chr><int><int><chr>
\n" ], "text/latex": [ - "A data.table: 6 × 4\n", + "A data.table: 6 \u00d7 4\n", "\\begin{tabular}{llll}\n", " \\#chr & start & end & gene\\_id\\\\\n", " & & & \\\\\n", @@ -590,7 +590,7 @@ ], "text/markdown": [ "\n", - "A data.table: 6 × 4\n", + "A data.table: 6 \u00d7 4\n", "\n", "| #chr <chr> | start <int> | end <int> | gene_id <chr> |\n", "|---|---|---|---|\n", @@ -722,7 +722,7 @@ " --genotype-file output/genotype_by_chrom/protocol_example.genotype.merged.plink_qc.genotype_by_chrom_files.txt \\\n", " --phenotype-file output/phenotype/phenotype_by_chrom_for_cis/bulk_rnaseq.phenotype_by_chrom_files.txt \\\n", " --covariate-file output/covariate/protocol_example.rnaseq.bed.protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.Marchenko_PC.gz \\\n", - " --cwd output/tensorqtl_cis --name protocol_example --MAC 5 --numThreads 2" + " --cwd output/tensorqtl_cis --cohort protocol_example --context bulk_rnaseq --modality mRNA --MAC 5 --numThreads 2" ] }, { @@ -737,12 +737,12 @@ "\n", " **Computational strategy designed for trans analysis:**\n", " \n", - " Trans analysis faces significant memory challenges as we calculate all associations between all molecular traits × all genetic variants across the genome, creating a massive computational burden. To address this challenge, we implement a two-stage chromosome-based parallelization approach:\n", + " Trans analysis faces significant memory challenges as we calculate all associations between all molecular traits \u00d7 all genetic variants across the genome, creating a massive computational burden. To address this challenge, we implement a two-stage chromosome-based parallelization approach:\n", "\n", " **Stage 1 (trans_1): Chromosome-based parallelization**\n", " - Phenotype data is processed per chromosome (e.g., 22 separate jobs for autosomes)\n", " - For each phenotype chromosome, we test associations against variants from all 22 chromosomes\n", - " - This creates phenotype_chr × genotype_chr combinations (e.g., phenotype chr1 vs genotype chr1-22); Garbage was collected between each chromosome combination caculation to release memory\n", + " - This creates phenotype_chr \u00d7 genotype_chr combinations (e.g., phenotype chr1 vs genotype chr1-22); Garbage was collected between each chromosome combination caculation to release memory\n", " - Results are combined across all chromosome combinations and saved as compressed files\n", "\n", " **Stage 2 (trans_2): Significance filtering**\n", @@ -773,7 +773,7 @@ " --genotype-file output/genotype_by_chrom/protocol_example.genotype.merged.plink_qc.genotype_by_chrom_files.txt \\\n", " --phenotype-file output/phenotype/phenotype_by_chrom_for_cis/bulk_rnaseq.phenotype_by_chrom_files.txt \\\n", " --covariate-file output/covariate/protocol_example.rnaseq.bed.protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.Marchenko_PC.gz \\\n", - " --cwd output/tensorqtl_trans --name protocol_example --MAC 5 --numThreads 2 \\\n", + " --cwd output/tensorqtl_trans --cohort protocol_example --context bulk_rnaseq --modality mRNA --MAC 5 --numThreads 2 \\\n", " --trans-geno-chromosome 22 --region-list data/combined_AD_genes.csv --region-list-phenotype-column 4" ] }, @@ -933,7 +933,10 @@ "# Optional pattern to filter covariates (list of covariate prefixes or exact names)\n", "parameter: covariate_pattern = []\n", "# Prefix for the analysis output\n", - "parameter: name = \"\"\n", + "parameter: cohort = str\n", + "parameter: context = str\n", + "parameter: modality = str\n", + "name = f\"{cohort}.{context}.{modality}\"\n", "# An optional subset of regions of molecular features to analyze. The last column is the gene names\n", "parameter: region_list = path()\n", "parameter: region_list_phenotype_column = 4\n", diff --git a/code/SoS/data_preprocessing/genotype/GWAS_QC.ipynb b/code/SoS/data_preprocessing/genotype/GWAS_QC.ipynb index 1710a8d7f..c2e20a1e9 100644 --- a/code/SoS/data_preprocessing/genotype/GWAS_QC.ipynb +++ b/code/SoS/data_preprocessing/genotype/GWAS_QC.ipynb @@ -149,6 +149,7 @@ "sos run pipeline/GWAS_QC.ipynb qc_no_prune \\\n", " --cwd output/gwas_qc/plink \\\n", " --genoFile output/genotype_formatting/plink/protocol_example.genotype.merged.bed \\\n", + " --cohort protocol_example \\\n", " --geno-filter 0.1 \\\n", " --mind-filter 0.1 \\\n", " --hwe-filter 1e-08 \\\n", @@ -187,7 +188,8 @@ "sos run pipeline/GWAS_QC.ipynb genotype_phenotype_sample_overlap \\\n", " --cwd output/gwas_qc/genotype \\\n", " --genoFile output/gwas_qc/plink/protocol_example.genotype.merged.plink_qc.fam \\\n", - " --phenoFile input/rnaseq/protocol_example.rnaseq.bed.gz\n" + " --phenoFile input/rnaseq/protocol_example.rnaseq.bed.gz \\\n", + " --cohort protocol_example\n" ] }, { @@ -224,7 +226,7 @@ "sos run pipeline/GWAS_QC.ipynb king \\\n", " --cwd output/gwas_qc/kinship \\\n", " --genoFile output/gwas_qc/plink/protocol_example.genotype.merged.plink_qc.bed \\\n", - " --name protocol_example.king \\\n", + " --cohort protocol_example.king \\\n", " --keep-samples output/gwas_qc/genotype/protocol_example.rnaseq.bed.sample_genotypes.txt\n" ] }, @@ -262,6 +264,7 @@ "sos run pipeline/GWAS_QC.ipynb qc \\\n", " --cwd output/gwas_qc/genotype \\\n", " --genoFile output/gwas_qc/kinship/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.bed \\\n", + " --cohort protocol_example \\\n", " --mac-filter 5\n" ] }, @@ -296,6 +299,7 @@ "sos run pipeline/GWAS_QC.ipynb qc \\\n", " --cwd output/gwas_qc/genotype \\\n", " --genoFile output/gwas_qc/plink/protocol_example.genotype.merged.plink_qc.bed \\\n", + " --cohort protocol_example \\\n", " --mac-filter 5\n" ] }, @@ -362,7 +366,8 @@ "# the output directory for generated files\n", "parameter: cwd = path(\"output\")\n", "# A string to identify your analysis run\n", - "parameter: name = \"\"\n", + "parameter: cohort = str\n", + "name = cohort\n", "# PLINK binary files (either BED/BIM/FAM or PGEN/PVAR/PSAM format)\n", "parameter: genoFile = paths\n", "# The path to the file that contains the list of samples to remove (format FID, IID)\n", diff --git a/code/SoS/data_preprocessing/genotype/PCA.ipynb b/code/SoS/data_preprocessing/genotype/PCA.ipynb index 19178d94b..e27163ee9 100644 --- a/code/SoS/data_preprocessing/genotype/PCA.ipynb +++ b/code/SoS/data_preprocessing/genotype/PCA.ipynb @@ -80,7 +80,7 @@ "source": [ "### **Step 1.** Estimate kinship in the sample (prerequisite)\n", "\n", - "Before PCA we need to know which individuals are related, because PCA must be computed on an **unrelated** subset and the related individuals are projected back afterwards. This step runs the `king` workflow from GWAS_QC to estimate pairwise kinship and split the samples into *unrelated* and *related* sets. It is an upstream prerequisite — its outputs (the `king.unrelated` and `king.related` PLINK bundles) feed the PCA steps below." + "Before PCA we need to know which individuals are related, because PCA must be computed on an **unrelated** subset and the related individuals are projected back afterwards. This step runs the `king` workflow from GWAS_QC to estimate pairwise kinship and split the samples into *unrelated* and *related* sets. It is an upstream prerequisite \u2014 its outputs (the `king.unrelated` and `king.related` PLINK bundles) feed the PCA steps below." ] }, { @@ -210,7 +210,7 @@ "kernel": "SoS" }, "source": [ - "**What you should see:** the workflow runs `flashpca_1` (PCA), then `detect_outliers` (Mahalanobis distance per population) and `plot_pca`, finishing with \"executed successfully with 5 completed steps\". The PCA model and PC scores are written to the `--cwd` directory. The scree plot below shows the proportion of variance explained by each PC — use it to decide how many PCs to keep." + "**What you should see:** the workflow runs `flashpca_1` (PCA), then `detect_outliers` (Mahalanobis distance per population) and `plot_pca`, finishing with \"executed successfully with 5 completed steps\". The PCA model and PC scores are written to the `--cwd` directory. The scree plot below shows the proportion of variance explained by each PC \u2014 use it to decide how many PCs to keep." ] }, { @@ -232,7 +232,7 @@ "source": [ "### **Step 4.** Project the related individuals onto the PCA model\n", "\n", - "The related individuals were held out of the PCA, so we now project them onto the model fitted from the unrelated samples. First reduce the related genotypes to **exactly** the same pruned variant set used to build the model (the `plink2 --extract` cell below) — otherwise the projection errors with \"Input number of variants should be the same as used in the previous PCA model\". Then run `project_samples`, passing a phenotype table (must contain an `IID` column) and, if available, a population/ethnicity column via `--label-col` / `--pop-col` (here the `race` column) so the projected samples can be coloured and the outlier detection can run." + "The related individuals were held out of the PCA, so we now project them onto the model fitted from the unrelated samples. First reduce the related genotypes to **exactly** the same pruned variant set used to build the model (the `plink2 --extract` cell below) \u2014 otherwise the projection errors with \"Input number of variants should be the same as used in the previous PCA model\". Then run `project_samples`, passing a phenotype table (must contain an `IID` column) and, if available, a population/ethnicity column via `--label-col` / `--pop-col` (here the `race` column) so the projected samples can be coloured and the outlier detection can run." ] }, { @@ -253,6 +253,7 @@ " --pca-model output/pca_uf/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.plink_qc.prune.protocol_example.pca.rds \\\n", " --label-col race \\\n", " --pop-col race \\\n", + " --cohort protocol_example \\\n", " --maha-k 2" ] }, @@ -338,11 +339,11 @@ "kernel": "SoS" }, "source": [ - "## Optional — Analysis by population (admixed / multi-ancestry data)\n", + "## Optional \u2014 Analysis by population (admixed / multi-ancestry data)\n", "\n", - "The Steps 1–4 above produce a single PCA for the whole cohort, which is appropriate for a **homogeneous** population. If your cohort contains **multiple ancestry groups**, run the per-population workflow below instead: split the samples by population, then repeat the same QC → PCA → projection steps **separately within each population**.\n", + "The Steps 1\u20134 above produce a single PCA for the whole cohort, which is appropriate for a **homogeneous** population. If your cohort contains **multiple ancestry groups**, run the per-population workflow below instead: split the samples by population, then repeat the same QC \u2192 PCA \u2192 projection steps **separately within each population**.\n", "\n", - "The sub-steps below mirror Steps 1–4 exactly, just looped over each population (`for i in race1 race3`). Follow them in order:\n", + "The sub-steps below mirror Steps 1\u20134 exactly, just looped over each population (`for i in race1 race3`). Follow them in order:\n", "\n", "| Sub-step | What it does |\n", "| --- | --- |\n", @@ -350,7 +351,7 @@ "| **P1. QC unrelated** | Variant/sample QC within each population (unrelated samples) |\n", "| **P2. Extract variants in related** | Keep the P1 pruned variants in the related samples |\n", "| **P3. PCA on unrelated** | Run `flashpca` per population |\n", - "| **P4. Project related** | Project related samples onto each population’s PCA |\n", + "| **P4. Project related** | Project related samples onto each population\u2019s PCA |\n", "| **P5. Finalize** | Per-population QC to finalize genotypes |" ] }, @@ -386,9 +387,9 @@ "kernel": "SoS" }, "source": [ - "### P1 — QC on unrelated samples, per population\n", + "### P1 \u2014 QC on unrelated samples, per population\n", "\n", - "Same as Step 2 (QC), but looped over each population. Below we only show Populations 1 and 3 as an example. **Next:** once each population has a QC’d unrelated set, extract those variants from the related samples (P2)." + "Same as Step 2 (QC), but looped over each population. Below we only show Populations 1 and 3 as an example. **Next:** once each population has a QC\u2019d unrelated set, extract those variants from the related samples (P2)." ] }, { @@ -419,11 +420,11 @@ "kernel": "SoS" }, "source": [ - "### P2 — Extract the selected variants from related samples, per population\n", + "### P2 \u2014 Extract the selected variants from related samples, per population\n", "\n", "Same as the extract part of Step 4, looped per population (sample-level missingness filter only).\n", "\n", - "> Note: Population 1 may have **no related samples**, in which case its `*.related.for_pca_race1.filtered.extracted.stderr` log will say `Error: No people remaining after --keep.` That is expected — just continue with the populations that do have related samples. **Next:** run PCA on each population’s unrelated set (P3)." + "> Note: Population 1 may have **no related samples**, in which case its `*.related.for_pca_race1.filtered.extracted.stderr` log will say `Error: No people remaining after --keep.` That is expected \u2014 just continue with the populations that do have related samples. **Next:** run PCA on each population\u2019s unrelated set (P3)." ] }, { @@ -454,9 +455,9 @@ "kernel": "SoS" }, "source": [ - "### P3 — Run PCA (flashpca) on unrelated samples, per population\n", + "### P3 \u2014 Run PCA (flashpca) on unrelated samples, per population\n", "\n", - "Same as Step 3, looped per population. Here the number of PCs is set to 5 (`--k 5`) because each population’s sample size is small. **Next:** project the related samples onto each population’s PCA model (P4)." + "Same as Step 3, looped per population. Here the number of PCs is set to 5 (`--k 5`) because each population\u2019s sample size is small. **Next:** project the related samples onto each population\u2019s PCA model (P4)." ] }, { @@ -472,7 +473,7 @@ "source": [ "for i in 3; do\n", " sos run pipeline/PCA.ipynb flashpca \\\n", - " --name pop_$i \\\n", + " --cohort pop_$i \\\n", " --cwd output/pca_uf \\\n", " --genoFile output/pca_uf/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.plink_qc.pop_$i.plink_qc.prune.bed \\\n", " --phenoFile input/covariate/protocol_example.pca_pheno.txt \\\n", @@ -489,7 +490,7 @@ "kernel": "SoS" }, "source": [ - "### P4 — Project related samples onto each population’s PCA, per population\n", + "### P4 \u2014 Project related samples onto each population\u2019s PCA, per population\n", "\n", "Same as Step 4 (projection), looped per population. Only populations that actually have related samples are projected (here, Population 3 only). **Next:** finalize genotypes per population (P5)." ] @@ -507,7 +508,7 @@ "source": [ "for i in 3; do\n", " sos run pipeline/PCA.ipynb project_samples \\\n", - " --name pop_$i \\\n", + " --cohort pop_$i \\\n", " --cwd output/pca_uf \\\n", " --genoFile output/pca_uf/protocol_example.genotype.merged.plink_qc.protocol_example.king.related.for_pca.plink_qc.extracted.pop_$i.plink_qc.extracted.bed \\\n", " --phenoFile input/covariate/protocol_example.pca_pheno.txt \\\n", @@ -546,7 +547,7 @@ "kernel": "SoS" }, "source": [ - "### P5 — Finalize genotypes per population\n", + "### P5 \u2014 Finalize genotypes per population\n", "\n", "Same as the **Finalize genotype QC by PCA** step for a homogeneous population, applied per population and taking the detected outliers into account. See the `GWAS_QC.ipynb` documentation for the available QC options and recommendations. This is the end of the per-population workflow." ] @@ -584,7 +585,8 @@ "# the output directory for generated files\n", "parameter: cwd = path(\"output\")\n", "# A string to identify your analysis run\n", - "parameter: name = \"\"\n", + "parameter: cohort = str\n", + "name = cohort\n", "# Name of the population column in the phenoFile\n", "parameter: pop_col = \"\"\n", "# Name of the populations (from the population column) you would like to plot and show on the PCA plot\n", @@ -687,7 +689,7 @@ " data_rows = rows[1:] if has_header else rows\n", " if not data_rows or not all(len(r) >= 2 and r[0] == r[1] for r in data_rows):\n", " return path(pheno_path)\n", - " # FID==IID everywhere — write a normalized copy to out_dir with FID=0\n", + " # FID==IID everywhere \u2014 write a normalized copy to out_dir with FID=0\n", " out_path = os.path.join(str(out_dir), os.path.basename(pheno_path))\n", " with open(out_path, 'w', newline='') as fh:\n", " w = csv.writer(fh, delimiter='\\t')\n", diff --git a/code/SoS/data_preprocessing/genotype/genotype_formatting.ipynb b/code/SoS/data_preprocessing/genotype/genotype_formatting.ipynb index e2290602c..d16bf4785 100644 --- a/code/SoS/data_preprocessing/genotype/genotype_formatting.ipynb +++ b/code/SoS/data_preprocessing/genotype/genotype_formatting.ipynb @@ -118,6 +118,7 @@ "sos run pipeline/genotype_formatting.ipynb vcf_to_plink \\\n", " --genoFile `ls input/genotype/protocol_example.genotype.chr*.vcf.gz | grep -vE \"rawchr|withfmt|add_chr\"` \\\n", " --cwd output/genotype_formatting/plink \\\n", + " --cohort protocol_example \\\n", " -j 4\n" ] }, @@ -145,7 +146,7 @@ "kernel": "SoS" }, "source": [ - "This step merges all per-chromosome PLINK files into one genome-wide bundle. Note: only `plink` (v1.9) can perform the merge — `plink2` does not support it.\n" + "This step merges all per-chromosome PLINK files into one genome-wide bundle. Note: only `plink` (v1.9) can perform the merge \u2014 `plink2` does not support it.\n" ] }, { @@ -248,7 +249,8 @@ "parameter: numThreads = 20\n", "# the path to a bed file or VCF file, a vector of bed files or VCF files, or a text file listing the bed files or VCF files to process\n", "parameter: genoFile = paths\n", - "parameter: name = \"\"\n", + "parameter: cohort = str\n", + "name = cohort\n", "# use this function to edit memory string for PLINK input\n", "from sos.utils import expand_size\n", "cwd = f\"{cwd:a}\"\n", diff --git a/code/SoS/misc/genotype_alignment.ipynb b/code/SoS/misc/genotype_alignment.ipynb index f939fb025..8182c2709 100644 --- a/code/SoS/misc/genotype_alignment.ipynb +++ b/code/SoS/misc/genotype_alignment.ipynb @@ -91,7 +91,7 @@ "sos run pipeline/genotype_alignment.ipynb genotype_alignment \\\n", " --geno_list_paths input/genotype/protocol_example.geno_cohortA input/genotype/protocol_example.geno_cohortB \\\n", " --cwd output/genotype_alignment \\\n", - " --name protocol_example" + " --cohort protocol_example" ] }, { @@ -147,7 +147,8 @@ "import pandas as pd\n", "## Path to work directory where output locates\n", "parameter: cwd = path(\"output\")\n", - "parameter: name = \"demo\"\n", + "parameter: cohort = str\n", + "name = cohort\n", "## Containers that contains the necessary packages\n", "parameter: container = \"\"\n", "# For cluster jobs, number commands to run per job\n", @@ -286,4 +287,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} \ No newline at end of file +} From a8bbd14b9ac0d4fc0ff86488562460aef0818fd0 Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 8 Jul 2026 03:32:51 -0400 Subject: [PATCH 02/12] fix(enrichment): correct eoo_enrichment CLI params and sldsc postprocessing function name eoo_enrichment's documented example passed --name, which the [global] step never declared (it only had parameter: name = 'eoo'); replaced with the real required parameters --trait/--annotation-name so the example runs as written. sldsc_enrichment's docs/text referenced pecotmr::sldsc_postprocessing_pipeline, which does not exist in the installed pecotmr version; the real exported function is sldscPostprocessingPipeline. Updated all references. Note: the postprocess/meta_subset steps still fail on a deeper argument-signature mismatch (old flat-argument calling convention vs. the new S4-object-based SldscData API) which is not fixed here and needs a follow-up. --- code/SoS/enrichment/eoo_enrichment.ipynb | 9 ++++++--- code/SoS/enrichment/sldsc_enrichment.ipynb | 14 +++++++------- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/code/SoS/enrichment/eoo_enrichment.ipynb b/code/SoS/enrichment/eoo_enrichment.ipynb index 343ddec3b..14de6e95a 100644 --- a/code/SoS/enrichment/eoo_enrichment.ipynb +++ b/code/SoS/enrichment/eoo_enrichment.ipynb @@ -303,7 +303,8 @@ "sos run pipeline/eoo_enrichment.ipynb enrichment \\\n", " --significant_variants_path input/enrichment/protocol_example.eoo_significant_variants.tsv.gz \\\n", " --baseline_anno_path input/enrichment/protocol_example.eoo_baseline_annotation.tsv \\\n", - " --name protocol_example \\\n", + " --trait protocol_example \\\n", + " --annotation-name baseline \\\n", " --cwd output/eoo_enrichment" ] }, @@ -366,7 +367,9 @@ "# Number of threads\n", "parameter: numThreads = 8\n", "# For cluster jobs, number commands to run per job\n", - "parameter: name = 'eoo'\n", + "parameter: trait = str\n", + "parameter: annotation_name = str\n", + "name = f\"{trait}.{annotation_name}\"\n", "parameter: job_size = 1\n", "parameter: walltime = '12h'\n", "parameter: mem = '16G'" @@ -612,4 +615,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} \ No newline at end of file +} diff --git a/code/SoS/enrichment/sldsc_enrichment.ipynb b/code/SoS/enrichment/sldsc_enrichment.ipynb index 8b352789c..9821381df 100644 --- a/code/SoS/enrichment/sldsc_enrichment.ipynb +++ b/code/SoS/enrichment/sldsc_enrichment.ipynb @@ -10,7 +10,7 @@ "\n", "Minimal working-example driver for the S-LDSC functional-enrichment pipeline. The **Steps** section below gives one ready-to-run `sos run` command per workflow, using the toy inputs symlinked under `input/`.\n", "\n", - "> **Environment note.** Steps 1–2 (`make_annotation_files_ldscore`, `get_heritability`) wrap the external **polyfun** toolkit (`compute_ldscores.py`, `ldsc.py`, `munge_polyfun_sumstats.py`) and require pre-computed reference-panel files (baseline-LD scores, LD weights, `.frq`, and PLINK `.bed/.bim/.fam`). polyfun is **not installed in this environment** and the reference panel is not shipped with the toy example, so those two steps cannot be executed here; their commands are provided for use on a system where polyfun and a matching panel are available. Steps 3–4 (`postprocess`, `meta_subset`) use `pecotmr::sldsc_postprocessing_pipeline` (available here) and read the `.results`/`.log` files produced by Step 2.\n" + "> **Environment note.** Steps 1–2 (`make_annotation_files_ldscore`, `get_heritability`) wrap the external **polyfun** toolkit (`compute_ldscores.py`, `ldsc.py`, `munge_polyfun_sumstats.py`) and require pre-computed reference-panel files (baseline-LD scores, LD weights, `.frq`, and PLINK `.bed/.bim/.fam`). polyfun is **not installed in this environment** and the reference panel is not shipped with the toy example, so those two steps cannot be executed here; their commands are provided for use on a system where polyfun and a matching panel are available. Steps 3–4 (`postprocess`, `meta_subset`) use `pecotmr::sldscPostprocessingPipeline` (available here) and read the `.results`/`.log` files produced by Step 2.\n" ] }, { @@ -48,7 +48,7 @@ "\n", "**Stage 1 - polyfun.** Three SoS workflows wrap polyfun: `make_annotation_files_ldscore` converts target annotations into polyfun `.annot.gz` and runs `compute_ldscores.py` (toggles `compute_single` and `compute_joint`, both default `True`; the joint dir is only emitted when $N \\geq 2$); `munge_sumstats_polyfun` preprocesses each GWAS into LDSC format; `get_heritability` runs polyfun's `ldsc.py` once per `--target-anno-dir`, enforcing the MAF cutoff via `--frqfile-chr` (`maf_cutoff` accepts only `0` or `0.05`).\n", "\n", - "**Stage 2 - pecotmr post-processing.** A single `pecotmr::sldsc_postprocessing_pipeline` call consumes all polyfun outputs: it extracts $\\tau$, $E$, $h^2_g$, EnrichStat p-value and per-block jackknife $\\tau$ values; computes $sd_C$ and $M_{\\mathrm{ref}}$ over the regression's MAF-cutoff SNP set; standardizes $\\tau \\to \\tau^*$ for single and joint modes; auto-detects binary vs continuous annotations; and runs a DerSimonian-Laird random-effects meta-analysis across traits, producing three meta tables ($\\tau^*$ cross-type comparable, $E$ within-binary, EnrichStat within-binary). Output is an R list with `per_trait` and `meta` entries.\n", + "**Stage 2 - pecotmr post-processing.** A single `pecotmr::sldscPostprocessingPipeline` call consumes all polyfun outputs: it extracts $\\tau$, $E$, $h^2_g$, EnrichStat p-value and per-block jackknife $\\tau$ values; computes $sd_C$ and $M_{\\mathrm{ref}}$ over the regression's MAF-cutoff SNP set; standardizes $\\tau \\to \\tau^*$ for single and joint modes; auto-detects binary vs continuous annotations; and runs a DerSimonian-Laird random-effects meta-analysis across traits, producing three meta tables ($\\tau^*$ cross-type comparable, $E$ within-binary, EnrichStat within-binary). Output is an R list with `per_trait` and `meta` entries.\n", "\n", "**Stage 3 - subset meta-analysis.** `pecotmr::meta_sldsc_random` re-runs the meta on a trait subset without re-running the regression (lightweight, interactive):\n", "\n", @@ -180,7 +180,7 @@ "- EnrichStat point estimate and its standard error (formula below) — **(pecotmr: `standardize_sldsc_trait`)**.\n", "- DerSimonian-Laird random-effects meta-analysis of $\\tau^*_C$, $E_C$, or EnrichStat across traits — **(pecotmr: `meta_sldsc_random`)**.\n", "\n", - "The top-level entry point `pecotmr::sldsc_postprocessing_pipeline` orchestrates all of the above.\n", + "The top-level entry point `pecotmr::sldscPostprocessingPipeline` orchestrates all of the above.\n", "\n", "#### Standardized tau ($\\tau^*$) — (pecotmr)\n", "\n", @@ -484,7 +484,7 @@ "source": [ "## Step 3. `Post-processing (pecotmr) and meta-analysis`\n", "\n", - "*Post-Processing (`pecotmr::sldsc_postprocessing_pipeline`)*\n", + "*Post-Processing (`pecotmr::sldscPostprocessingPipeline`)*\n", "\n", "A single R function call consumes all polyfun outputs for the run and produces the final tables:\n", "\n", @@ -499,7 +499,7 @@ "\n", "The `[postprocess]` step reads all polyfun outputs under `heritability_cwd`\n", "(which contains the $N$ single-target subdirectories and optionally the\n", - "joint subdirectory) and calls `pecotmr::sldsc_postprocessing_pipeline()`\n", + "joint subdirectory) and calls `pecotmr::sldscPostprocessingPipeline()`\n", "to produce per-trait standardized tables and the default random-effects\n", "meta across all traits.\n", "\n", @@ -1304,7 +1304,7 @@ "outputs": [], "source": [ "[postprocess]\n", - "# Post-processing of polyfun outputs via pecotmr::sldsc_postprocessing_pipeline.\n", + "# Post-processing of polyfun outputs via pecotmr::sldscPostprocessingPipeline.\n", "# Reads .results / .log / .part_delete for all traits in `traits_file`, both\n", "# single-target and (when present) joint-target runs, computes Gazal-style\n", "# tau*, EnrichStat with back-solved jackknife SE, and runs the default\n", @@ -1357,7 +1357,7 @@ " trait_joint_prefix <- setNames(rep(NA_character_, length(traits)), traits)\n", " }\n", "\n", - " res <- sldsc_postprocessing_pipeline(\n", + " res <- sldscPostprocessingPipeline(\n", " trait_single_prefixes = trait_single_prefixes,\n", " trait_joint_prefix = trait_joint_prefix,\n", " target_anno_dir = \"${target_anno_dir}\",\n", From 5db31458359e7bb4b62128dec4130bb704e851e4 Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 8 Jul 2026 03:33:18 -0400 Subject: [PATCH 03/12] fix(xqtl_modifier_score): fix config paths, silent-failure guard, and prediction model-path default - Default --config-dir pointed at code/xqtl_modifier_score, which does not exist; the real location is code/SoS/xqtl_modifier_score. - The [train] bash step ran 'touch $[_output]' unconditionally after the python call with no 'set -e', so a failed training run was still reported as successful. Added set -e to the block. - data_config.yaml used ../.. for paths relative to the script, but the script actually runs from code/SoS/xqtl_modifier_score (3 levels deep, not 2), so every input path resolved incorrectly. Fixed to ../../... - ems_prediction's documented --model-path pointed at output/ems_training/protocol_example_chr2_scEEMS_model.joblib, which is never produced; corrected to the real path under output/xqtl_modifier_score/protocol_example/model_results/. Validated end-to-end: train produces a real cross-validated model (joblib), predict produces real per-variant EMS scores. --- code/SoS/xqtl_modifier_score/data_config.yaml | 14 +++++++------- .../xqtl_modifier_score/ems_prediction.ipynb | 12 ++++++------ code/SoS/xqtl_modifier_score/ems_training.ipynb | 17 +++++++++-------- 3 files changed, 22 insertions(+), 21 deletions(-) diff --git a/code/SoS/xqtl_modifier_score/data_config.yaml b/code/SoS/xqtl_modifier_score/data_config.yaml index a79e01bd8..b3b1bfda7 100644 --- a/code/SoS/xqtl_modifier_score/data_config.yaml +++ b/code/SoS/xqtl_modifier_score/data_config.yaml @@ -2,7 +2,7 @@ feature_data: gene_constraint: name: "gene_lof" - file_path: "../../input/xqtl_modifier_score/protocol_example.gene_constraint.xlsx" + file_path: "../../../input/xqtl_modifier_score/protocol_example.gene_constraint.xlsx" xlsx_sheet: "Supplementary Table 1" include_columns: ["ensg", "post_mean"] column_mapping: @@ -15,22 +15,22 @@ feature_data: population_genetics: name: "maf" - file_pattern: "../../input/xqtl_modifier_score/protocol_example.gnomad_MAF_chr{chromosome}.tsv" + file_pattern: "../../../input/xqtl_modifier_score/protocol_example.gnomad_MAF_chr{chromosome}.tsv" column_mapping: variant_id: "variant_id" target_value: "gnomad_MAF" distance_features: - columns_dict_file: "../../input/xqtl_modifier_score/protocol_example.columns_dict.pkl" + columns_dict_file: "../../../input/xqtl_modifier_score/protocol_example.columns_dict.pkl" subset_keys: ["distance"] columns_to_remove: ["abs_distance_TSS", "distance_TSS"] regulatory_features: - columns_dict_file: "../../input/xqtl_modifier_score/protocol_example.columns_dict.pkl" + columns_dict_file: "../../../input/xqtl_modifier_score/protocol_example.columns_dict.pkl" subset_keys: ["ABC", "celltype", "baseline"] deep_learning_features: - columns_dict_file: "../../input/xqtl_modifier_score/protocol_example.columns_dict.pkl" + columns_dict_file: "../../../input/xqtl_modifier_score/protocol_example.columns_dict.pkl" subset_keys: ["chrombpnet_positive", "diff", "tf_positive"] transformations: absolute_value: ["diff", "tf_positive", "chrombpnet_positive"] @@ -39,13 +39,13 @@ feature_data: generated_columns: ["length_diff", "is_SNP", "is_indel", "is_insertion", "is_deletion", "gene_lof", "gnomad_MAF"] training_data: - base_dir: "../../input/xqtl_modifier_score/{cohort}/training_data" + base_dir: "../../../input/xqtl_modifier_score/{cohort}/training_data" file_pattern: "annotated_data_{cohort}_{chromosome}.parquet" train_dir_pattern: "train_NPR_{npr_tr}_PIP_{pos_threshold}_{neg_threshold}" test_dir_pattern: "test_NPR_{npr_te}_PIP_{pos_threshold}_{neg_threshold}" metadata_columns: ["variant_id", "pip", "CHR", "BP", "REF", "ALT", "SNP", "label", "weight"] output: - base_dir: "../../output/xqtl_modifier_score/{cohort}/model_results" + base_dir: "../../../output/xqtl_modifier_score/{cohort}/model_results" predictions_dir: "predictions_parquet_catboost" diff --git a/code/SoS/xqtl_modifier_score/ems_prediction.ipynb b/code/SoS/xqtl_modifier_score/ems_prediction.ipynb index 26f9ac82f..23b9c9a5e 100644 --- a/code/SoS/xqtl_modifier_score/ems_prediction.ipynb +++ b/code/SoS/xqtl_modifier_score/ems_prediction.ipynb @@ -34,7 +34,7 @@ "### Step 2: Execute Prediction Pipeline\n", "\n", "```bash\n", - "cd ~/xqtl-protocol/code/xqtl_modifier_score/\n", + "cd ~/xqtl-protocol/code/SoS/xqtl_modifier_score/\n", "python model_training_model5_only.py Mic_mega_eQTL 2 \\\n", " --data_config data_config.yaml \\\n", " --model_config model_config.yaml\n", @@ -400,8 +400,8 @@ "sos run pipeline/ems_prediction.ipynb predict \\\n", " --cohort protocol_example \\\n", " --chromosome 2 \\\n", - " --model_path output/ems_training/protocol_example_chr2_scEEMS_model.joblib \\\n", - " --data_config code/xqtl_modifier_score/data_config.yaml \\\n", + " --model_path output/xqtl_modifier_score/protocol_example/model_results/model_standard_subset_weighted_chr_chr2_NPR_1.joblib \\\n", + " --data_config code/SoS/xqtl_modifier_score/data_config.yaml \\\n", " --cwd output/ems_prediction\n" ] }, @@ -453,11 +453,11 @@ "# Chromosome to score\n", "parameter: chromosome = '2'\n", "# Trained CatBoost model (.joblib) from the EMS Training workflow\n", - "parameter: model_path = path('output/ems_training/protocol_example_chr2_scEEMS_model.joblib')\n", + "parameter: model_path = path('output/xqtl_modifier_score/protocol_example/model_results/model_standard_subset_weighted_chr_chr2_NPR_1.joblib')\n", "# Data configuration YAML (cohort, variant list, feature set)\n", - "parameter: data_config = path('code/xqtl_modifier_score/data_config.yaml')\n", + "parameter: data_config = path('code/SoS/xqtl_modifier_score/data_config.yaml')\n", "# Directory holding gems_pipeline.py\n", - "parameter: pipeline_dir = path('code/xqtl_modifier_score')\n", + "parameter: pipeline_dir = path('code/SoS/xqtl_modifier_score')\n", "parameter: job_size = 1\n", "parameter: mem = '60G'\n", "parameter: walltime = '24h'\n" diff --git a/code/SoS/xqtl_modifier_score/ems_training.ipynb b/code/SoS/xqtl_modifier_score/ems_training.ipynb index 6b687532a..81dc39fe8 100644 --- a/code/SoS/xqtl_modifier_score/ems_training.ipynb +++ b/code/SoS/xqtl_modifier_score/ems_training.ipynb @@ -189,7 +189,7 @@ "\n", "### Step 1: Running the GEMS Pipeline\n", "```bash\n", - "cd ~/xqtl-protocol/code/xqtl_modifier_score/\n", + "cd ~/xqtl-protocol/code/SoS/xqtl_modifier_score/\n", "python gems_pipeline.py train Mic_mega_eQTL 2 \\\n", " --data_config data_config.yaml \\\n", " --model_config model_config.yaml\n", @@ -216,8 +216,8 @@ "sos run pipeline/ems_training.ipynb train \\\n", " --cohort protocol_example \\\n", " --chromosome 2 \\\n", - " --data-config code/xqtl_modifier_score/data_config.yaml \\\n", - " --model-config code/xqtl_modifier_score/model_config.yaml \\\n", + " --data-config code/SoS/xqtl_modifier_score/data_config.yaml \\\n", + " --model-config code/SoS/xqtl_modifier_score/model_config.yaml \\\n", " --cwd output/ems_training\n" ] }, @@ -326,7 +326,7 @@ }, "outputs": [], "source": [ - "sos run code/xqtl_modifier_score/ems_training.ipynb -h" + "sos run code/SoS/xqtl_modifier_score/ems_training.ipynb -h" ] }, { @@ -357,11 +357,11 @@ "# Chromosome to train on\n", "parameter: chromosome = '2'\n", "# Data configuration YAML (cohort, eQTL paths, feature set)\n", - "parameter: data_config = path('code/xqtl_modifier_score/data_config.yaml')\n", + "parameter: data_config = path('code/SoS/xqtl_modifier_score/data_config.yaml')\n", "# Model configuration YAML (CatBoost hyperparameters, feature weighting)\n", - "parameter: model_config = path('code/xqtl_modifier_score/model_config.yaml')\n", + "parameter: model_config = path('code/SoS/xqtl_modifier_score/model_config.yaml')\n", "# Directory holding gems_pipeline.py\n", - "parameter: pipeline_dir = path('code/xqtl_modifier_score')\n", + "parameter: pipeline_dir = path('code/SoS/xqtl_modifier_score')\n", "parameter: job_size = 1\n", "parameter: mem = '60G'\n", "parameter: walltime = '24h'" @@ -379,6 +379,7 @@ "output: f'{cwd:a}/{cohort}_chr{chromosome}_scEEMS_model.done'\n", "task: trunk_workers = 1, trunk_size = job_size, mem = mem, walltime = walltime, tags = f'{step_name}_{cohort}_chr{chromosome}'\n", "bash: expand = \"$[ ]\", workdir = pipeline_dir, stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout'\n", + " set -e\n", " python gems_pipeline.py train $[cohort] $[chromosome] \\\n", " --data_config $[data_config:a] \\\n", " --model_config $[model_config:a]\n", @@ -447,7 +448,7 @@ "source": [ "### Using Your Trained Models\n", "\n", - "Once training is complete, load the trained models for predictions. Please refer to **[EMS Predictions](https://statfungen.github.io/xqtl-protocol/code/xqtl_modifier_score/ems_prediction.html)** for detailed prediction workflows and variant scoring." + "Once training is complete, load the trained models for predictions. Please refer to **[EMS Predictions](https://statfungen.github.io/xqtl-protocol/code/SoS/xqtl_modifier_score/ems_prediction.html)** for detailed prediction workflows and variant scoring." ] }, { From 5b13245e9592a231df04d5b48d58a175e806bf6b Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 8 Jul 2026 03:33:36 -0400 Subject: [PATCH 04/12] fix(reference_data): remove trailing dot in rss_ld_sketch cohort-id; fix ld_prune_reference variable name rss_ld_sketch's documented --cohort-id protocol_example. had a stray trailing dot, producing double-dot output filenames (protocol_example..pgen etc). Removed the trailing dot from both occurrences and regenerated the chr22 outputs to confirm clean single-dot filenames. ld_prune_reference referenced an undefined variable 'genotype' instead of the actual parameter 'genotype_list' when reading the genotype file list. --- code/SoS/reference_data/ld_prune_reference.ipynb | 2 +- code/SoS/reference_data/rss_ld_sketch.ipynb | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/code/SoS/reference_data/ld_prune_reference.ipynb b/code/SoS/reference_data/ld_prune_reference.ipynb index 544fbacb3..d16d77e89 100644 --- a/code/SoS/reference_data/ld_prune_reference.ipynb +++ b/code/SoS/reference_data/ld_prune_reference.ipynb @@ -160,7 +160,7 @@ "parameter: r2 = 1e-3\n", "\n", "import pandas as pd\n", - "geno_path = pd.read_csv(genotype, sep = \"\\t\")\n", + "geno_path = pd.read_csv(genotype_list, sep = \"\\t\")\n", "input_df = geno_path.values.tolist()\n", "input_blocks = [x[0] for x in input_df]\n", "input_files = [x[1:] for x in input_df] \n", diff --git a/code/SoS/reference_data/rss_ld_sketch.ipynb b/code/SoS/reference_data/rss_ld_sketch.ipynb index 0396c03d0..2e3e3ded4 100644 --- a/code/SoS/reference_data/rss_ld_sketch.ipynb +++ b/code/SoS/reference_data/rss_ld_sketch.ipynb @@ -158,7 +158,7 @@ " --output-dir output/rss_ld_sketch \\\n", " --W-matrix output/rss_ld_sketch/W_B50.npy \\\n", " --B 50 \\\n", - " --cohort-id protocol_example. \\\n", + " --cohort-id protocol_example \\\n", " --cwd output/rss_ld_sketch" ] }, @@ -240,7 +240,7 @@ "source": [ "sos run pipeline/rss_ld_sketch.ipynb merge_chrom \\\n", " --output-dir output/rss_ld_sketch \\\n", - " --cohort-id protocol_example. \\\n", + " --cohort-id protocol_example \\\n", " --chrom 22 \\\n", " --cwd output/rss_ld_sketch" ] From f7c846af3ddaaa1bf96daf31a77831eaf626a89c Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 8 Jul 2026 03:33:54 -0400 Subject: [PATCH 05/12] fix(mnm-and-pecotmr): parameter/API fixes across TWAS and colocalization notebooks - mnm_regression: drop unavailable susieInf method from default twas_methods. - mnm_postprocessing: rename unused 'name' placeholder to required --study, actually wired into the output name. - univariate_fine_mapping_twas_vignette: documented command skipped the required qtl_dataset_construct step and pointed at a phenotype file that does not match the multi-trait manifest schema; fixed to chain qtl_dataset_construct+susie_twas with the correct manifest-schema files. - SuSiE_enloc: fixed analysis_name_prefix matching, which used split('.')[0] assuming the wrong filename format; corrected to a case-insensitive substring match against the real {prefix}.{cohort}.{analysis_name}.{gene}...rds naming; also corrected the documented enrichment output structure/example to match real output. - twas_ctwas: added a fallback that writes explicit empty placeholder outputs when a region produces none, so downstream steps do not break on missing files. - legacy_twas_weights_convert.R: fixed output field name cvPerformance -> cvResult to match the consumer's expected schema. --- .../mnm_methods/mnm_regression.ipynb | 2 +- .../SoS/mnm_analysis/mnm_postprocessing.ipynb | 9 ++++---- ...nivariate_fine_mapping_twas_vignette.ipynb | 20 ++++++++++------- .../SoS/pecotmr_integration/SuSiE_enloc.ipynb | 22 ++++++++++++++++--- code/SoS/pecotmr_integration/twas_ctwas.ipynb | 14 +++++++----- .../legacy_twas_weights_convert.R | 2 +- 6 files changed, 46 insertions(+), 23 deletions(-) diff --git a/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb b/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb index 075cd8a8a..308326297 100644 --- a/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb +++ b/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb @@ -711,7 +711,7 @@ "# legacy regional_data + susie_twas.R path.\n", "parameter: cis_window = 1000000\n", "parameter: fine_mapping_methods = \"susie\"\n", - "parameter: twas_methods = \"susie,susieInf,mrash,enet,lasso,mcp,scad,l0learn,bayes_r,bayes_c\"\n", + "parameter: twas_methods = \"susie,mrash,enet,lasso,mcp,scad,l0learn,bayes_r,bayes_c\"\n", "parameter: fine_mapping_coverage = 0.95\n", "# Comma-separated context names to restrict both passes to; empty = all contexts.\n", "parameter: contexts = \"\"\n", diff --git a/code/SoS/mnm_analysis/mnm_postprocessing.ipynb b/code/SoS/mnm_analysis/mnm_postprocessing.ipynb index 44848a570..79929fc3f 100644 --- a/code/SoS/mnm_analysis/mnm_postprocessing.ipynb +++ b/code/SoS/mnm_analysis/mnm_postprocessing.ipynb @@ -80,7 +80,7 @@ "# {file_path}/{prefix}.{region}.{suffix}\n", "# example: ROSMAP_AC_eQTL.ENSG00000012779.univariate_bvsr.rds\n", "# \\__ prefix __/ \\___ region ___/ \\___ suffix ___/\n", - "# --name tag used in output filenames and intermediate RDS names\n", + "# --study tag used in output filenames and intermediate RDS names\n", "# (final output: {cwd}/summary/{name}.{qtl_type}.top_loci.bed.gz)\n", "# --qtl_type cis | trans | trans_ appears in the combined filename\n", "# --region_file bed of regions (chr, start, end, region_id) to enumerate\n", @@ -97,7 +97,7 @@ " --file_path /path/to/fine_mapping \\\n", " --prefix \\\n", " --suffix \\\n", - " --name \\\n", + " --study \\\n", " --min_corr 0.8 \\\n", " --geno_ref ./Fungen_xQTL.ROSMAP_NIA_WGS.ROSMAP_NIA_WGS.MSBB_WGS_ADSP_hg38.MiGA.MAP_Brain-xQTL_Gwas_geno_0.STARNET.aligned.bim.gz \\\n", " --qtl_type cis \\\n", @@ -130,7 +130,7 @@ "outputs": [], "source": [ "sos run pipeline/mnm_postprocessing.ipynb susie_to_tsv \\\n", - " --cwd output_postproc --name toy \\\n", + " --cwd output_postproc --study toy \\\n", " --rds_path input_postproc/cs.*.uni.bvsr.rds \\\n", " --region-list input_postproc/region_list.tsv \\\n", " --container \"\" -j1" @@ -211,7 +211,8 @@ "import pandas as pd\n", "# A region list file documenting the chr_pos_ref_alt of a susie_object\n", "parameter: cwd = path(\"output\")\n", - "parameter: name = \"demo\"\n", + "parameter: study = str\n", + "name = study\n", "\n", "## Path to work directory where output locates\n", "## Containers that contains the necessary packages\n", diff --git a/code/SoS/mnm_analysis/univariate_fine_mapping_twas_vignette.ipynb b/code/SoS/mnm_analysis/univariate_fine_mapping_twas_vignette.ipynb index ca8cfc3e4..d1e8ef26f 100644 --- a/code/SoS/mnm_analysis/univariate_fine_mapping_twas_vignette.ipynb +++ b/code/SoS/mnm_analysis/univariate_fine_mapping_twas_vignette.ipynb @@ -127,14 +127,18 @@ "outputs": [], "source": [ "# Peak/chromatin example (same pipeline, different input)\n", - "sos run pipeline/mnm_regression.ipynb susie_twas \\\n", - " --name protocol_example \\\n", - " --cwd output_10peaks \\\n", - " --genoFile input/genotype/protocol_example.genotype.chr22.bed \\\n", - " --phenoFile input/phenotype/protocol_example.pheno_manifest.tsv \\\n", - " --covFile input/covariate/protocol_example.covariates.tsv \\\n", - " --customized-association-windows input/finemapping/protocol_example.association_windows.bed \\\n", - " -j1\n" + "# NOTE: this is a peak catalog with one phenotype file per peak; the\n", + "# qtl_dataset_construct manifest schema (ID, path, cond) instead expects a\n", + "# single multi-trait phenotype file per context. protocol_example.peaks_merged.bed.gz\n", + "# (all 10 peaks, built from peaks_split/) and protocol_example.pheno_manifest_context.tsv\n", + "# match that schema; run qtl_dataset_construct before susie_twas (as in Step 1 above).\n", + "sos run pipeline/mnm_regression.ipynb qtl_dataset_construct+susie_twas \\\n", + "--name protocol_example --cwd output/susie_twas_peaks \\\n", + "--genoFile input/genotype/protocol_example.genotype.chr22.bed \\\n", + "--phenoFile input/phenotype/protocol_example.pheno_manifest_context.tsv \\\n", + "--covFile input/covariate/protocol_example.covariates.tsv \\\n", + "--customized-association-windows input/finemapping/protocol_example.association_windows.bed \\\n", + "--region-name C22P107555 -j1\n" ] }, { diff --git a/code/SoS/pecotmr_integration/SuSiE_enloc.ipynb b/code/SoS/pecotmr_integration/SuSiE_enloc.ipynb index c560d6fd8..d10d51339 100644 --- a/code/SoS/pecotmr_integration/SuSiE_enloc.ipynb +++ b/code/SoS/pecotmr_integration/SuSiE_enloc.ipynb @@ -124,7 +124,23 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Output\n\nOutput file: `{cwd}/{name}.{context}.enrichment.rds`\n\n```r\nenrich <- readRDS(\"output/output/xqtl_gwas_enrichment/protocol_example.enloc.Knight_eQTL_brain.enrichment.rds\")\n\nlength(enrich) # 2\nnames(enrich) # [1] \"\" \"unused_xqtl_variants\"\n\nstr(enrich[[1]])\n# List of 10\n# $ Alternative (coloc) p1 : num 0.000558 prior prob. a random SNP is a causal xQTL variant\n# $ Alternative (coloc) p12 : num 6.31e-08 prior prob. a random SNP is causal for both xQTL and GWAS\n# $ Alternative (coloc) p2 : num 0.000113 prior prob. a random SNP is a causal GWAS variant\n# $ Effective MI rounds : num 25 number of multiple-imputation iterations\n# $ Enrichment (no shrinkage): num -0.241 enrichment slope a1 (raw)\n# $ Enrichment (w/ shrinkage): num -0.000102 enrichment slope a1 (shrinkage-regularised, used for p12)\n# $ Intercept : num -7.49 baseline log-odds a0\n# $ sd (intercept) : num 0.316 SE of a0\n# $ sd (no shrinkage) : num 48.6 SE of a1 (raw)\n# $ sd (w/ shrinkage) : num 1 SE of a1 (shrinkage)\n\nlength(enrich$unused_xqtl_variants[[1]]) # 18\nhead(enrich$unused_xqtl_variants[[1]])\n# [1] \"1:20568238:C:T\" \"1:20621145:G:A\" \"1:21427833:C:G\" ...\n```\n\n> **Note on toy data**: estimates are produced but unreliable (1 gene \u00d7 1 GWAS region). `sd (no shrinkage)` = 48.6 reflects near-zero information; shrinkage regularises `a1` to \u2248 0. Reliable a0/a1 require thousands of gene \u00d7 GWAS region pairs." + "### Output\n", + "Output file: `{cwd}/{name}.{context}.enrichment.rds`\n", + "\n", + "```r\n", + "enrich <- readRDS(\"output/xqtl_gwas_enrichment/protocol_example.enloc.protocol_example.enloc.Knight_eQTL_brain.enrichment.rds\")\n", + "\n", + "str(enrich)\n", + "# data.frame: 1 obs. of 6 variables:\n", + "# $ gwasStudy : chr \"AD_Bellenguez_2022\"\n", + "# $ qtlStudy : chr \"protocol_example.enloc.protocol_example.enloc\"\n", + "# $ qtlContext : chr \"Knight_eQTL_brain\"\n", + "# $ enrichment : num NA\n", + "# $ enrichmentSe : num NA\n", + "# $ enrichmentLogOdds : num NA\n", + "```\n", + "\n", + "> **Note on toy data**: one row per gwas/qtl-context pair. Enrichment values are `NA` here because a single gene x single GWAS block does not give the underlying logistic-regression enrichment model enough data to fit; this is expected for the toy MWE, not a pipeline error. Reliable enrichment estimates require many gene x GWAS-region pairs." ] }, { @@ -323,7 +339,7 @@ " \"\"\"Merge and filter the dataframes based on condition/context where context is a substring of condition, and analysis_name to pick the corresponding original files.\"\"\"\n", " # Explode the QTL data into separate rows\n", " new_df = new_df.set_index(['condition', 'region_id', 'GWAS_original_data'])['QTL_original_data'].str.split(',', expand=True).stack().reset_index(name='QTL_original_data').drop('level_3', axis=1)\n", - " new_df['analysis_name_prefix'] = new_df['QTL_original_data'].apply(lambda x: x.split('.')[0])\n", + " new_df['analysis_name_prefix'] = new_df['QTL_original_data'] # bugfix: keep full filename; matched via substring below, not position[0], since real file names are '{prefix}.{cohort}.{analysis_name}.{gene}...rds' not '{analysis_name}...rds'\n", "\n", " # Create a custom merge logic to match context as a substring of condition\n", " def custom_merge(row, df_meta):\n", @@ -339,7 +355,7 @@ " merged_df = pd.concat(results, keys=new_df.index).reset_index(level=1, drop=True).join(new_df, how='outer')\n", "\n", " # Filter rows where analysis_name_prefix matches analysis_name exactly\n", - " filtered_df = merged_df[merged_df['analysis_name_prefix'].str.strip() == merged_df['analysis_name'].str.strip()]\n", + " filtered_df = merged_df[merged_df.apply(lambda r: str(r['analysis_name']).strip().lower() in str(r['analysis_name_prefix']).lower(), axis=1)] # bugfix: case-insensitive substring match instead of exact-equality on a mis-parsed prefix\n", "\n", " # Drop unnecessary columns and duplicates\n", " filtered_df = filtered_df.drop(columns=['analysis_name_prefix', 'context']).drop_duplicates()\n", diff --git a/code/SoS/pecotmr_integration/twas_ctwas.ipynb b/code/SoS/pecotmr_integration/twas_ctwas.ipynb index 4640fefc6..10acc761b 100644 --- a/code/SoS/pecotmr_integration/twas_ctwas.ipynb +++ b/code/SoS/pecotmr_integration/twas_ctwas.ipynb @@ -33,7 +33,7 @@ "\n", "*Notes:* specifying a `column_file_path` (YAML) enables `load_rss_data()`, which standardizes column names and can generate missing `z`/`beta` columns via `col_to_flip`; without it, the simpler `tabix_region()` is used. Use `--comment_string \"#\"` to handle comment lines in the summary statistics file (by default no comment symbol is assumed). For LASSO, Elastic Net, and mr.ash the weights are taken as-is for QTL variants overlapping GWAS variants, while SuSiE weights can be adjusted to exactly match GWAS variants.\n", "\n", - "**Input — GWAS data interface (similar to `susie_rss`):**\n", + "**Input \u2014 GWAS data interface (similar to `susie_rss`):**\n", "- GWAS summary statistics files: tab-delimited, tabix-indexed by `chrom`/`pos`; first 4 columns `chrom`, `pos`, `A1`, `A2`. For MR, `effect_allele_frequency` and sample-size columns are required.\n", "- Optional column-mapping YAML: required `chrom`, `pos`, `A1`, `A2`, `z` (or `betahat`/`sebetahat`); optional `n`, `var_y`.\n", "- Optional GWAS meta-file (`study_id`, chrom, file path, optional mapping file; chrom `0` = genome-wide):\n", @@ -49,7 +49,7 @@ "chr1 0 6480000 ENSG00000008128 1724356 KNIGHT_pQTL.ENSG00000008128.univariate_susie_twas_weights.rds, ... (truncated)\n", "```\n", "\n", - "**Output — TWAS result table** (imputable genes only): `gwas_study, chrom, start, end, block, gene, TSS, context, is_imputable, method, is_selected_method, rsq_adj_cv, pval_cv, twas_z, twas_pval`. Key columns: `TSS` = transcription start site; `start`/`end` = gene window from the [extended TADB window](https://github.com/cumc/xqtl-analysis/blob/main/resource/TADB_enhanced_cis.coding.bed); `is_imputable` = CV r-square >0.01 and pvalue <0.05 in >=1 method; `is_selected_method` = best model (highest CV r-square, CV pvalue <0.05); `block` = LD region of the gene's TSS based on [predefined LD blocks](https://github.com/cumc/xqtl-data/blob/main/descriptor/reference_data/ld_reference.md).\n", + "**Output \u2014 TWAS result table** (imputable genes only): `gwas_study, chrom, start, end, block, gene, TSS, context, is_imputable, method, is_selected_method, rsq_adj_cv, pval_cv, twas_z, twas_pval`. Key columns: `TSS` = transcription start site; `start`/`end` = gene window from the [extended TADB window](https://github.com/cumc/xqtl-analysis/blob/main/resource/TADB_enhanced_cis.coding.bed); `is_imputable` = CV r-square >0.01 and pvalue <0.05 in >=1 method; `is_selected_method` = best model (highest CV r-square, CV pvalue <0.05); `block` = LD region of the gene's TSS based on [predefined LD blocks](https://github.com/cumc/xqtl-data/blob/main/descriptor/reference_data/ld_reference.md).\n", "```\n", "chr molecular_id TSS start end context gwas_study method is_imputable is_selected_method rsq_cv pval_cv twas_z twas_pval\n", "1 ENSG00000060718 103108871 101000000 104000000 AC_DeJager_eQTL Bellenguez_EADB_2022 bayes_r TRUE TRUE 0.25 3.39e-39 0.39 0.69\n", @@ -172,13 +172,13 @@ "kernel": "SoS" }, "source": [ - "## Step 2: cTWAS — region assembly, global parameters, and fine-mapping\n", + "## Step 2: cTWAS \u2014 region assembly, global parameters, and fine-mapping\n", "\n", "Combine the selected best-performing TWAS prediction models, TWAS gene z-scores, GWAS SNP z-scores, and LD reference region information into the per-region cTWAS input; estimate the global group prior across all regions; then screen and fine-map regions to obtain posterior inclusion probabilities (PIP) for genes and SNPs. The section is split into three runs: region-data assembly, global-parameter estimation (`--run_param_est`), and causal fine-mapping (`--run_finemapping`).\n", "\n", "The cTWAS R package is required: `remotes::install_github(\"xinhe-lab/ctwas\", ref = \"multigroup\")`. Variant selection is performed based on the `top_loci` table, SuSiE CS set, or twas-weights effect size when any of `twas_weight_cutoff` (default=0), `cs_min_cor` (default=0), `min_pip_cutoff` (default=0), or `max_num_variants` (default=Inf) is set away from its default. The prior can use a single shared group via `--prior_var_structure shared_all` or the multi-group model via `--multi_group`.\n", "\n", - "**Input — TWAS region information** (multiple TWAS and SNP data within each region are combined for joint inference to select variables — genes or SNPs — likely to be directly associated with the phenotype rather than via correlation):\n", + "**Input \u2014 TWAS region information** (multiple TWAS and SNP data within each region are combined for joint inference to select variables \u2014 genes or SNPs \u2014 likely to be directly associated with the phenotype rather than via correlation):\n", "```\n", "chrom start end block_id\n", "1 1000 5000 block1\n", @@ -186,7 +186,7 @@ "3 3000 7000 block3\n", "```\n", "\n", - "**Output — cTWAS fine-mapping** (per variable: `id, molecular_id, type, context, susie_pip, group, cs, region_id, z`):\n", + "**Output \u2014 cTWAS fine-mapping** (per variable: `id, molecular_id, type, context, susie_pip, group, cs, region_id, z`):\n", "```\n", "id molecular_id type context susie_pip cs region_id z\n", "ENSG00000148429|eQTL_Inh_DeJager_eQTL ENSG00000148429 eQTL Inh_DeJager_eQTL 0.07121779052573 L1 10_10500888_12817813 -6.20856901762341\n", @@ -920,6 +920,7 @@ "\n", " if (length(weight_db_list_update) == 0 || all(sapply(weight_db_list_update, length) == 0)) {\n", " message(\"No valid twas weight files found after filtering. Exiting the script.\")\n", + " for (of in c(${_output:r,})) { if (grepl(\"[.]rds$\", of)) saveRDS(list(), of) else readr::write_tsv(data.frame(message = character(0)), of) }\n", " quit(save = \"no\", status = 0) \n", " }\n", "\n", @@ -967,6 +968,7 @@ " # Merging with xQTL meta-data \n", " if (is.null(twas_results_db$twas_result) || nrow(twas_results_db$twas_result) == 0) {\n", " message(\"twas_results_db$twas_result is NULL. Exiting script normally.\")\n", + " for (of in c(${_output:r,})) { if (grepl(\"[.]rds$\", of)) saveRDS(list(), of) else readr::write_tsv(data.frame(message = character(0)), of) }\n", " quit(save = \"no\", status = 0) \n", " }\n", "\n", @@ -1036,4 +1038,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file diff --git a/code/script/pecotmr_integration/legacy_twas_weights_convert.R b/code/script/pecotmr_integration/legacy_twas_weights_convert.R index 1050621d6..a17792da1 100644 --- a/code/script/pecotmr_integration/legacy_twas_weights_convert.R +++ b/code/script/pecotmr_integration/legacy_twas_weights_convert.R @@ -67,7 +67,7 @@ for (trait in names(legacy)) { wvec <- stats::setNames(as.numeric(wmat), vids) cv <- perfToCv(perfL[[paste0(tok, "_performance")]]) entries[[length(entries) + 1L]] <- TwasWeightsEntry( - variantIds = vids, weights = wvec, cvPerformance = cv, + variantIds = vids, weights = wvec, cvResult = cv, standardized = FALSE, dataType = context) rs <- c(rs, argv$study); rc <- c(rc, context) rt <- c(rt, trait); rm <- c(rm, tok) From b5b01d54f9a89bfeeeec6245de9cb2353287947d Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 8 Jul 2026 03:34:11 -0400 Subject: [PATCH 06/12] fix(molecular-phenotypes): support gene-ID peaks without positional info in pseudobulk BED construction to_midpoint_bed assumed every peak ID encodes chr-start-end and could be split on '-', which fails for gene-ID-based peaks (e.g. plain Ensembl IDs). Added a GTF-based gene-position lookup (load_gene_positions) used as a fallback to build the BED record (TSS-based single-base coordinate) when the peak ID isn't itself a coordinate triplet. snRNAseq_preprocessing: no functional change, only re-saved with escaped unicode in comments (box-drawing separators / arrows). --- .../QC/pseudobulk_preprocessing.ipynb | 158 +++++++++++------- .../snRNAseq_preprocessing.ipynb | 50 +++--- 2 files changed, 124 insertions(+), 84 deletions(-) diff --git a/code/SoS/molecular_phenotypes/QC/pseudobulk_preprocessing.ipynb b/code/SoS/molecular_phenotypes/QC/pseudobulk_preprocessing.ipynb index 6ee25a83c..05be36280 100644 --- a/code/SoS/molecular_phenotypes/QC/pseudobulk_preprocessing.ipynb +++ b/code/SoS/molecular_phenotypes/QC/pseudobulk_preprocessing.ipynb @@ -72,14 +72,14 @@ "\n", "### Process\n", "\n", - "1. Load each Seurat object and subset to target cell type — skips objects where cell type is not present\n", + "1. Load each Seurat object and subset to target cell type \u2014 skips objects where cell type is not present\n", "2. Merge all subsets across objects and join layers\n", "3. Aggregate raw counts by sample (`AggregateExpression`)\n", "4. Filter out samples with fewer than `min_cells` cells (default: 10)\n", - "5. Strip Ensembl version suffixes from gene IDs (`ENSG00000000010.1` → `ENSG00000000010`)\n", - "6. Save as `pseudobulk_counts_{celltype}.csv.gz` — raw counts only, normalization handled downstream in `pseudobulk_qc`\n", + "5. Strip Ensembl version suffixes from gene IDs (`ENSG00000000010.1` \u2192 `ENSG00000000010`)\n", + "6. Save as `pseudobulk_counts_{celltype}.csv.gz` \u2014 raw counts only, normalization handled downstream in `pseudobulk_qc`\n", "\n", - "> **GLU cell type**: Due to its large size, process in two batches (files 1–6 and 7–11) and pass separately.\n", + "> **GLU cell type**: Due to its large size, process in two batches (files 1\u20136 and 7\u201311) and pass separately.\n", "\n", "### Parameters\n", "\n", @@ -94,10 +94,10 @@ "\n", "| File | Description |\n", "|------|-------------|\n", - "| `pseudobulk_counts_{celltype}.csv.gz` | Raw pseudobulk count matrix (genes × samples) |\n", + "| `pseudobulk_counts_{celltype}.csv.gz` | Raw pseudobulk count matrix (genes \u00d7 samples) |\n", "\n", "\n", - "**Timing:** 10–30 min per cell type depending on object size" + "**Timing:** 10\u201330 min per cell type depending on object size" ] }, { @@ -129,27 +129,27 @@ "\n", "| File | Description |\n", "|------|-------------|\n", - "| `rosmap_sample_mapping_data.csv` | Mapping reference: `individualID → sampleid` |\n", + "| `rosmap_sample_mapping_data.csv` | Mapping reference: `individualID \u2192 sampleid` |\n", "| `metadata_{celltype}.csv` | Per-cell-type sample metadata |\n", "| `pseudobulk_peaks_counts_{celltype}.csv.gz` *(snATAC-seq)* | Per-cell-type peak count matrices |\n", "| `pseudobulk_counts_{celltype}.csv.gz` *(snRNA-seq)* | Per-cell-type gene count matrices |\n", "\n", "### Process\n", "\n", - "**Part 1 — Metadata files**\n", + "**Part 1 \u2014 Metadata files**\n", "\n", "For each metadata file:\n", "1. Look up each `individualID` in the mapping reference\n", - "2. Assign `sampleid` — falls back to `individualID` if no mapping found\n", + "2. Assign `sampleid` \u2014 falls back to `individualID` if no mapping found\n", "3. Reorder columns: `sampleid` first, then `individualID`, then the rest\n", "4. Save updated file\n", "\n", - "**Part 2 — Count matrix files**\n", + "**Part 2 \u2014 Count matrix files**\n", "\n", "For each count file:\n", "1. Extract the header row (column names only)\n", "2. Keep the first column (peak or gene IDs) unchanged\n", - "3. Map remaining column names (`individualID` → `sampleid`) where mapping exists, otherwise keep original\n", + "3. Map remaining column names (`individualID` \u2192 `sampleid`) where mapping exists, otherwise keep original\n", "4. Write new header and stream data rows unchanged\n", "5. Recompress with gzip\n", "\n", @@ -157,7 +157,7 @@ "\n", "| Parameter | Default | Description |\n", "|-----------|---------|-------------|\n", - "| `map_file` | *required* | CSV with `individualID` → `sampleid` mapping |\n", + "| `map_file` | *required* | CSV with `individualID` \u2192 `sampleid` mapping |\n", "| `meta_files` | *required* | Metadata CSV files to remap |\n", "| `count_files` | *required* | Count CSV.gz files to remap |\n", "| `output_dir` | *required* | Parent output directory; writes to `{output_dir}/1_files_with_sampleid/` |\n", @@ -221,7 +221,7 @@ "8. Apply expression filtering (`filterByExpr`):\n", " - `min_count = 5`: minimum reads in at least one sample\n", " - `min_total_count = 15`: minimum total reads across all samples\n", - " - `min_prop = 0.1`: feature expressed in ≥10% of samples\n", + " - `min_prop = 0.1`: feature expressed in \u226510% of samples\n", "9. TMM normalization\n", "10. ***(Optional)*** Batch correction on `sequencingBatch`:\n", " - `limma::removeBatchEffect` (default)\n", @@ -238,7 +238,7 @@ "~ {tech_vars} + [sequencingBatch] + [Library]\n", "```\n", "> `sequencingBatch` and `Library` included only if present and have more than one level.\n", - "> Biological variables (`pmi`, `study`, `msex`, `age_death` etc.) are **not** included — they should not be regressed out as they may be associated with genotype.\n", + "> Biological variables (`pmi`, `study`, `msex`, `age_death` etc.) are **not** included \u2014 they should not be regressed out as they may be associated with genotype.\n", "\n", "### Parameters\n", "\n", @@ -350,7 +350,8 @@ "source": [ "sos run pipeline/pseudobulk_preprocessing.ipynb phenotype_formatting \\\n", " --residual-files input/snrnaseq/atac_residuals/MIC/protocol_example.snrnaseq.MIC_residuals.txt \\\n", - " --output-dir output/snrna_seq\n" + " --output-dir output/snrna_seq \\\n", + " --gtf-file input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf\n" ] }, { @@ -488,6 +489,7 @@ " Workflow Options:\n", " --residual-files (as list)\n", " --output-dir VAL (as str, required)\n", + " --gtf-file VAL (as path, required)\n", "```" ] }, @@ -532,7 +534,7 @@ "\n", " message(\"Loading and subsetting Seurat objects for: \", celltype)\n", "\n", - " # ── 1. Load and subset each Seurat object ─────────────────────\n", + " # \u2500\u2500 1. Load and subset each Seurat object \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " subsets <- list()\n", " for (f in seurat_files) {\n", " message(\"Loading: \", basename(f))\n", @@ -549,7 +551,7 @@ " if (length(subsets) == 0) stop(\"No Seurat objects contain celltype: \", celltype)\n", " message(\"Found \", celltype, \" in \", length(subsets), \" objects\")\n", "\n", - " # ── 2. Merge and aggregate ────────────────────────────────────\n", + " # \u2500\u2500 2. Merge and aggregate \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " merged <- Reduce(merge, subsets)\n", " rm(subsets)\n", " gc()\n", @@ -562,7 +564,7 @@ " rm(merged)\n", " gc()\n", "\n", - " # ── 3. Filter samples with < min_cells ────────────────────────\n", + " # \u2500\u2500 3. Filter samples with < min_cells \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " valid_samples <- names(cell_counts[cell_counts >= min_cells])\n", " # Seurat AggregateExpression replaces underscores in identities with dashes;\n", " # map valid_samples to the dash-form used in colnames(expr) for subsetting.\n", @@ -572,10 +574,10 @@ " colnames(expr) <- gsub(\"-\", \"_\", colnames(expr)) # restore underscore sample names\n", " message(\"Samples after min_cells (>= \", min_cells, \") filter: \", ncol(expr))\n", "\n", - " # ── 4. Strip Ensembl version suffixes ─────────────────────────\n", + " # \u2500\u2500 4. Strip Ensembl version suffixes \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " rownames(expr) <- gsub(\"\\\\..*$\", \"\", rownames(expr))\n", "\n", - " # ── 5. Save as csv.gz ─────────────────────────────────────────\n", + " # \u2500\u2500 5. Save as csv.gz \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " message(\"Genes: \", nrow(expr), \" | Samples: \", ncol(expr))\n", " dt <- data.table(gene_id = rownames(expr), as.data.frame(expr))\n", " fwrite(dt, out_file, compress = \"gzip\")\n", @@ -620,7 +622,10 @@ "import tempfile\n", "\n", "map_df = pd.read_csv(\"${map_file}\")\n", - "id_map = dict(zip(map_df[\"individualID\"], map_df[\"sampleid\"]))\n", + "if \"sampleid\" in map_df.columns:\n", + " id_map = dict(zip(map_df[\"individualID\"], map_df[\"sampleid\"]))\n", + "else:\n", + " id_map = dict(zip(map_df[\"individualID\"], map_df[\"individualID\"]))\n", "output_dir = \"${output_dir}/1_files_with_sampleid\"\n", "meta_files = ${meta_files}\n", "count_files = ${count_files}\n", @@ -642,7 +647,7 @@ " return str(val)\n", " return str(val)\n", "\n", - "# ── Process metadata ───────────────────────────────────────────────────────\n", + "# \u2500\u2500 Process metadata \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", "for in_path in meta_files:\n", " fname = os.path.basename(in_path)\n", " out_path = os.path.join(output_dir, fname)\n", @@ -661,7 +666,7 @@ " for _, row in meta.iterrows():\n", " writer.writerow([format_value(val) for val in row])\n", "\n", - "# ── Process count files ────────────────────────────────────────────────────\n", + "# \u2500\u2500 Process count files \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", "for in_path in count_files:\n", " fname = os.path.basename(in_path)\n", " out_path = os.path.join(output_dir, fname)\n", @@ -729,7 +734,7 @@ " library(GenomicRanges)\n", " if (as.logical(\"${batch_correction}\") && \"${batch_method}\" == \"combat\") library(sva)\n", "\n", - " # ── predictOffset ──────────────────────────────────────────────────────\n", + " # \u2500\u2500 predictOffset \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " predictOffset <- function(fit, tech_vars) {\n", " D <- fit$design\n", " Dm <- D\n", @@ -772,7 +777,7 @@ " parse_regions <- function(region_str) {\n", " if (is.null(region_str) || region_str == \"\") return(NULL)\n", " lapply(strsplit(region_str, \",\")[[1]], function(r) {\n", - " parts <- strsplit(trimws(r), \":|−|-\")[[1]]\n", + " parts <- strsplit(trimws(r), \":|\u2212|-\")[[1]]\n", " list(chr=parts[1], start=as.integer(parts[2]), end=as.integer(parts[3]))\n", " })\n", " }\n", @@ -798,7 +803,7 @@ " if (length(meta_files) != length(count_files))\n", " stop(\"meta_files and count_files must have the same length and order.\")\n", "\n", - " # ── Load tech vars from file ───────────────────────────────────────────\n", + " # \u2500\u2500 Load tech vars from file \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " tech_df <- fread(\"${tech_vars_file}\")\n", " tech_vars <- setdiff(colnames(tech_df), \"sampleid\")\n", " message(\"Tech vars: \", paste(tech_vars, collapse=\", \"))\n", @@ -821,19 +826,19 @@ " outdir <- file.path(\"${output_dir}/2_residuals\", ct)\n", " dir.create(outdir, recursive=TRUE, showWarnings=FALSE)\n", "\n", - " # ── 1. Load counts ─────────────────────────────────────────────────\n", + " # \u2500\u2500 1. Load counts \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " counts_raw <- fread(counts_file)\n", " counts <- as.matrix(counts_raw[, -1, with=FALSE])\n", " rownames(counts) <- counts_raw[[1]]\n", " rm(counts_raw)\n", "\n", - " # ── Auto-detect modality ───────────────────────────────────────────\n", + " # \u2500\u2500 Auto-detect modality \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " is_atac <- grepl(\"^chr.*-[0-9]+-[0-9]+$\", rownames(counts)[1])\n", " feat_label <- ifelse(is_atac, \"peaks\", \"genes\")\n", " message(\"Modality: \", ifelse(is_atac, \"snATAC-seq\", \"snRNA-seq\"))\n", " message(\"Loaded: \", nrow(counts), \" \", feat_label, \" x \", ncol(counts), \" samples\")\n", "\n", - " # ── 1b. Region/gene filtering (optional) ──────────────────────────\n", + " # \u2500\u2500 1b. Region/gene filtering (optional) \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " if (is_atac && !is.null(regions)) {\n", " message(\"Filtering peaks to specified regions...\")\n", " counts <- filter_regions(counts, regions)\n", @@ -844,25 +849,25 @@ " counts <- counts[genes_present, , drop=FALSE]\n", " }\n", "\n", - " # ── 2. Load metadata ───────────────────────────────────────────────\n", + " # \u2500\u2500 2. Load metadata \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " meta <- fread(meta_file)\n", " idcol <- intersect(c(\"sampleid\",\"sampleID\",\"individualID\",\"projid\"), colnames(meta))[1]\n", " if (is.na(idcol)) stop(\"Cannot find sample ID column in metadata.\")\n", "\n", - " # ── 3. Nuclei filter ──────────────────────────────────────────────\n", + " # \u2500\u2500 3. Nuclei filter \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " n_nuclei_col <- intersect(c(\"n_nuclei\",\"n.nuclei\",\"nNuclei\",\"nuclei_count\"), colnames(meta))[1]\n", " if (!is.na(n_nuclei_col)) {\n", " meta <- meta[meta[[n_nuclei_col]] > ${min_nuclei}]\n", " message(\"Samples after nuclei (>${min_nuclei}) filter: \", nrow(meta))\n", " }\n", "\n", - " # ── 4. Align samples ──────────────────────────────────────────────\n", + " # \u2500\u2500 4. Align samples \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " common <- intersect(meta[[idcol]], colnames(counts))\n", " if (length(common) == 0) stop(\"Zero sample overlap between metadata and count matrix.\")\n", " counts <- counts[, common, drop=FALSE]\n", " message(\"Samples after alignment: \", length(common))\n", "\n", - " # ── 5. Blacklist filtering ─────────────────────────────────────────\n", + " # \u2500\u2500 5. Blacklist filtering \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " if (\"${blacklist_file}\" != \"\" && file.exists(\"${blacklist_file}\")) {\n", " counts <- filter_blacklist(counts, \"${blacklist_file}\", feat_label)\n", " message(feat_label, \" after blacklist filter: \", nrow(counts))\n", @@ -870,17 +875,17 @@ " message(\"No blacklist file - skipping.\")\n", " }\n", "\n", - " # ── 6. Merge tech vars by sampleid ────────────────────────────────\n", + " # \u2500\u2500 6. Merge tech vars by sampleid \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " tech_sub <- tech_df[tech_df$sampleid %in% common]\n", " tech_sub <- tech_sub[match(common, tech_sub$sampleid)]\n", "\n", - " # ── 7. Drop samples with NA in tech vars ──────────────────────────\n", + " # \u2500\u2500 7. Drop samples with NA in tech vars \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " keep_rows <- complete.cases(tech_sub[, ..tech_vars])\n", " tech_sub <- tech_sub[keep_rows]\n", " counts <- counts[, tech_sub$sampleid, drop=FALSE]\n", " message(\"Valid samples for modelling: \", nrow(tech_sub))\n", "\n", - " # ── 8. Expression filtering ────────────────────────────────────────\n", + " # \u2500\u2500 8. Expression filtering \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " dge <- DGEList(counts=counts, samples=tech_sub)\n", " dge$samples$group <- factor(rep(\"all\", ncol(dge)))\n", " message(feat_label, \" before filter: \", nrow(dge))\n", @@ -896,10 +901,10 @@ " file.path(outdir, paste0(ct, \"_filtered_raw_counts.txt\")),\n", " sep=\"\\t\", quote=FALSE, col.names=NA)\n", "\n", - " # ── 9. TMM normalization ───────────────────────────────────────────\n", + " # \u2500\u2500 9. TMM normalization \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " dge <- calcNormFactors(dge, method=\"TMM\")\n", "\n", - " # ── 10. Optional batch correction ──────────────────────────────────\n", + " # \u2500\u2500 10. Optional batch correction \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " if (as.logical(\"${batch_correction}\") && \"sequencingBatch\" %in% colnames(dge$samples)) {\n", " batches <- dge$samples$sequencingBatch\n", " batch_counts <- table(batches)\n", @@ -922,7 +927,7 @@ " }\n", " }\n", "\n", - " # ── 11. Add batch vars to model if multi-level ────────────────────\n", + " # \u2500\u2500 11. Add batch vars to model if multi-level \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " batch_vars <- c()\n", " if (\"sequencingBatch\" %in% colnames(dge$samples) &&\n", " length(unique(dge$samples$sequencingBatch)) > 1) {\n", @@ -935,7 +940,7 @@ " batch_vars <- c(batch_vars, \"Library_factor\")\n", " }\n", "\n", - " # ── 12. Build design matrix ────────────────────────────────────────\n", + " # \u2500\u2500 12. Build design matrix \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " all_model_vars <- intersect(c(tech_vars, batch_vars), colnames(dge$samples))\n", " form <- as.formula(paste(\"~\", paste(all_model_vars, collapse=\" + \")))\n", " design <- model.matrix(form, data=dge$samples)\n", @@ -948,23 +953,23 @@ " }\n", " message(\"Design matrix: \", nrow(design), \" x \", ncol(design))\n", "\n", - " # ── 13. Voom + lmFit + eBayes ─────────────────────────────────────\n", + " # \u2500\u2500 13. Voom + lmFit + eBayes \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " v <- voom(dge, design, plot=FALSE)\n", " fit <- lmFit(v, design)\n", " fit <- eBayes(fit)\n", "\n", - " # ── 14. Offset + residuals ─────────────────────────────────────────\n", + " # \u2500\u2500 14. Offset + residuals \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " off <- predictOffset(fit, tech_vars=tech_vars)\n", " res <- residuals(fit, v$E)\n", " final <- off + res\n", "\n", - " # ── 15. Save residuals ─────────────────────────────────────────────\n", + " # \u2500\u2500 15. Save residuals \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " out_file <- file.path(outdir, paste0(ct, \"_residuals.txt\"))\n", " write.table(final, out_file, sep=\"\\t\", quote=FALSE, col.names=NA)\n", " message(\"Saved: \", out_file)\n", " message(\" \", ifelse(is_atac,\"Peaks\",\"Genes\"), \": \", nrow(final), \" | Samples: \", ncol(final))\n", "\n", - " # ── 16. Optional quantile normalization ───────────────────────────\n", + " # \u2500\u2500 16. Optional quantile normalization \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " if (as.logical(\"${quant_norm}\")) {\n", " final_qn <- t(apply(final, 1, rank, ties.method=\"average\"))\n", " final_qn <- stats::qnorm(final_qn / (ncol(final_qn) + 1))\n", @@ -1019,6 +1024,7 @@ "[phenotype_formatting]\n", "parameter: residual_files = []\n", "parameter: output_dir = str\n", + "parameter: gtf_file = path\n", "\n", "import os\n", "\n", @@ -1033,9 +1039,11 @@ " import os\n", " import subprocess\n", " import pandas as pd\n", + " import re\n", "\n", " residual_files = ${residual_files}\n", " output_dir = \"${output_dir}\"\n", + " gtf_file = \"${gtf_file}\"\n", "\n", " def read_residuals(path):\n", " first_line = open(path).readline().rstrip(\"\\n\")\n", @@ -1051,19 +1059,50 @@ " df.columns = col_names[1:]\n", " return peak_ids, df\n", "\n", - " def to_midpoint_bed(peak_ids, residuals):\n", - " parts = pd.Series(peak_ids).str.split(\"-\", expand=True)\n", - " chrs = parts[0].values\n", - " starts = parts[1].astype(int).values\n", - " ends = parts[2].astype(int).values\n", - " mids = ((starts + ends) // 2).astype(int)\n", - " bed = pd.DataFrame({\n", - " \"#chr\": chrs,\n", - " \"start\": mids,\n", - " \"end\": mids + 1,\n", - " \"ID\": peak_ids\n", - " })\n", - " bed = pd.concat([bed, residuals.reset_index(drop=True)], axis=1)\n", + " def load_gene_positions(gtf_path):\n", + " positions = {}\n", + " if not gtf_path:\n", + " return positions\n", + " with open(gtf_path) as fh:\n", + " for line in fh:\n", + " if line.startswith(\"#\"):\n", + " continue\n", + " fields = line.rstrip(\"\\n\").split(\"\\t\")\n", + " if len(fields) < 9 or fields[2] != \"gene\":\n", + " continue\n", + " chrom, start, end, strand, attrs = fields[0], int(fields[3]), int(fields[4]), fields[6], fields[8]\n", + " m = re.search(r'gene_id \"([^\"]+)\"', attrs)\n", + " if not m:\n", + " continue\n", + " tss = start if strand == \"+\" else end\n", + " positions[m.group(1)] = (chrom, tss)\n", + " return positions\n", + "\n", + " def to_midpoint_bed(peak_ids, residuals, gene_positions=None):\n", + " parts = pd.Series(peak_ids).str.split(\"-\", expand=True)\n", + " if parts.shape[1] >= 3:\n", + " chrs = parts[0].values\n", + " starts = parts[1].astype(int).values\n", + " ends = parts[2].astype(int).values\n", + " mids = ((starts + ends) // 2).astype(int)\n", + " bed = pd.DataFrame({\n", + " \"#chr\": chrs,\n", + " \"start\": mids,\n", + " \"end\": mids + 1,\n", + " \"ID\": peak_ids\n", + " })\n", + " bed = pd.concat([bed, residuals.reset_index(drop=True)], axis=1)\n", + " else:\n", + " gene_positions = gene_positions or {}\n", + " rows, keep = [], []\n", + " for idx, pid in enumerate(peak_ids):\n", + " gene_id = str(pid).split(\".\")[0]\n", + " if gene_id in gene_positions:\n", + " chrom, tss = gene_positions[gene_id]\n", + " rows.append((chrom, tss, tss + 1, pid))\n", + " keep.append(idx)\n", + " bed = pd.DataFrame(rows, columns=[\"#chr\", \"start\", \"end\", \"ID\"])\n", + " bed = pd.concat([bed, residuals.iloc[keep].reset_index(drop=True)], axis=1)\n", " return bed.sort_values([\"#chr\", \"start\"]).reset_index(drop=True)\n", "\n", " def run_cmd(cmd, label):\n", @@ -1075,6 +1114,7 @@ "\n", " out_dir = os.path.join(output_dir, \"3_pheno_reformat\")\n", " os.makedirs(out_dir, exist_ok=True)\n", + " gene_positions = load_gene_positions(gtf_file)\n", "\n", " for res_path in residual_files:\n", " ct = os.path.basename(os.path.dirname(res_path))\n", @@ -1088,7 +1128,7 @@ " peak_ids, residuals = read_residuals(res_path)\n", " print(f\"Loaded {len(peak_ids)} peaks x {residuals.shape[1]} samples\")\n", "\n", - " bed = to_midpoint_bed(peak_ids, residuals)\n", + " bed = to_midpoint_bed(peak_ids, residuals, gene_positions)\n", " out_bed = os.path.join(out_dir, f\"{ct}_phenotype.bed\")\n", " bed.to_csv(out_bed, sep=\"\\t\", index=False, float_format=\"%.15f\")\n", " print(f\"Written: {out_bed}\")\n", @@ -1128,4 +1168,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file diff --git a/code/SoS/molecular_phenotypes/snRNAseq_preprocessing.ipynb b/code/SoS/molecular_phenotypes/snRNAseq_preprocessing.ipynb index abe47cada..cfcb2650d 100644 --- a/code/SoS/molecular_phenotypes/snRNAseq_preprocessing.ipynb +++ b/code/SoS/molecular_phenotypes/snRNAseq_preprocessing.ipynb @@ -124,7 +124,7 @@ "import re\n", "from sos.utils import expand_size\n", "\n", - ") if container else \"\"\n", + "parameter: entrypoint= ('micromamba run -a \"\" -n' + ' ' + re.sub(r'(_apptainer:latest|_docker:latest|\\.sif)$', '', container.split('/')[-1])) if container else \"\"\n", "\n", "cwd = path(f'{cwd:a}')" ] @@ -187,11 +187,11 @@ " message(\"Doublet detection method : \", doublet_method)\n", " message(\"Ambient RNA removal method: \", ambient_method)\n", "\n", - " # ── Create output directories ─────────────────────────────────\n", + " # \u2500\u2500 Create output directories \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " dir.create(dirname(rds_out), showWarnings = FALSE, recursive = TRUE)\n", " dir.create(dirname(qc_table_out), showWarnings = FALSE, recursive = TRUE)\n", "\n", - " # ── QC summary table ──────────────────────────────────────────\n", + " # \u2500\u2500 QC summary table \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " qc_summary <- data.frame(Step = character(),\n", " Cells = numeric(), Genes = numeric(),\n", " stringsAsFactors = FALSE)\n", @@ -201,13 +201,13 @@ " Cells = ncol(obj), Genes = nrow(obj)))\n", " }\n", "\n", - " # ── Sample metadata ───────────────────────────────────────────\n", + " # \u2500\u2500 Sample metadata \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " samplelist <- read.csv(\"${sample_meta}\")\n", " samplelist$column_name <- paste0(samplelist$libraryBatch, \"-counts_\", samplelist$cellBarcode)\n", " samplelist <- samplelist[, c(\"individualID\", \"column_name\")]\n", " colnames(samplelist) <- c(\"sample\", \"column_name\")\n", "\n", - " # ── Helper: convert SCE to Seurat ─────────────────────────────\n", + " # \u2500\u2500 Helper: convert SCE to Seurat \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " convert_to_seurat <- function(sce) {\n", " obj <- CreateSeuratObject(counts = counts(sce),\n", " meta.data = as.data.frame(colData(sce)),\n", @@ -216,7 +216,7 @@ " return(obj)\n", " }\n", "\n", - " # ── Helper: assign sample IDs, convert back to SCE ────────────\n", + " # \u2500\u2500 Helper: assign sample IDs, convert back to SCE \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " process_sce <- function(seurat_obj, samplelist) {\n", " seurat_obj <- seurat_obj[, colnames(seurat_obj) %in% samplelist$column_name]\n", " update_summary(\"Removing Cells with no associated Patient ID\", seurat_obj)\n", @@ -229,16 +229,16 @@ " return(as.SingleCellExperiment(seurat_obj))\n", " }\n", "\n", - " # ── Helper: build QC algorithms list from chosen methods ──────\n", + " # \u2500\u2500 Helper: build QC algorithms list from chosen methods \u2500\u2500\u2500\u2500\u2500\u2500\n", " # Doublet detection methods (runCellQC YAML names):\n", - " # \"scds\" → cxds_bcds_hybrid (cxds + bcds hybrid score)\n", - " # \"scrublet\" → scrublet\n", - " # \"doubletFinder\" → doubletFinder\n", - " # \"doubletCells\" → doubletCells (scran-based)\n", - " # \"none\" → skip doublet detection\n", + " # \"scds\" \u2192 cxds_bcds_hybrid (cxds + bcds hybrid score)\n", + " # \"scrublet\" \u2192 scrublet\n", + " # \"doubletFinder\" \u2192 doubletFinder\n", + " # \"doubletCells\" \u2192 doubletCells (scran-based)\n", + " # \"none\" \u2192 skip doublet detection\n", " # Ambient RNA methods:\n", - " # \"decontX\" → decontX\n", - " # \"none\" → skip ambient RNA removal\n", + " # \"decontX\" \u2192 decontX\n", + " # \"none\" \u2192 skip ambient RNA removal\n", " build_algorithms <- function(doublet_method, ambient_method) {\n", " algos <- \"QCMetrics\"\n", " if (doublet_method == \"scds\") algos <- c(algos, \"cxds_bcds_hybrid\")\n", @@ -249,7 +249,7 @@ " return(algos)\n", " }\n", "\n", - " # ── Helper: build filter strings from chosen methods ──────────\n", + " # \u2500\u2500 Helper: build filter strings from chosen methods \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " build_filters <- function(doublet_method, ambient_method) {\n", " filters <- c(\n", " \"mito_percent < ${max_mito_percent}\",\n", @@ -271,7 +271,7 @@ " return(filters)\n", " }\n", "\n", - " # ── Helper: run SCTK-QC, filter, cluster ─────────────────────\n", + " # \u2500\u2500 Helper: run SCTK-QC, filter, cluster \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " process_sce_to_seurat <- function(sce) {\n", " algos <- build_algorithms(doublet_method, ambient_method)\n", " message(\"Running QC algorithms: \", paste(algos, collapse = \", \"))\n", @@ -283,20 +283,20 @@ " mitoGeneLocation = \"rownames\",\n", " seed = 12345)\n", "\n", - " # ── Write per-sample QC table ─────────────────────────────\n", + " # \u2500\u2500 Write per-sample QC table \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " sce_temp <- sampleSummaryStats(sce, sample = colData(sce)$sample, simple = FALSE)\n", " sst <- getSampleSummaryStatsTable(sce_temp, statsName = \"qc_table\")\n", " write.csv(sst, qc_table_out)\n", " rm(sce_temp)\n", "\n", - " # ── Apply filters ─────────────────────────────────────────\n", + " # \u2500\u2500 Apply filters \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " qc_filters <- build_filters(doublet_method, ambient_method)\n", " message(\"Applying filters:\")\n", " for (f in qc_filters) message(\" \", f)\n", " sce <- subsetSCECols(sce, colData = qc_filters)\n", " update_summary(\"SCTK-QC Results\", sce)\n", "\n", - " # ── Convert to Seurat and cluster ─────────────────────────\n", + " # \u2500\u2500 Convert to Seurat and cluster \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " seurat_obj <- CreateSeuratObject(counts = counts(sce),\n", " meta.data = as.data.frame(colData(sce)),\n", " min.cells = ncol(sce) * 0.005)\n", @@ -312,22 +312,22 @@ " return(seurat_obj)\n", " }\n", "\n", - " # ── 1. Import ─────────────────────────────────────────────────\n", + " # \u2500\u2500 1. Import \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " sce <- importCellRanger(cellRangerDirs = input_dir)\n", " update_summary(\"Raw Data\", sce)\n", "\n", - " # ── 2. Remove fake mapping genes ──────────────────────────────\n", + " # \u2500\u2500 2. Remove fake mapping genes \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " fake_genes <- grep(\"^(1|2|3|4|5|6|7|8|9|10|11|12|13|14|15|16|17|18|19|20|21|22|X|Y)_\",\n", " rownames(sce), value = TRUE)\n", " sce <- sce[setdiff(rownames(sce), fake_genes), ]\n", " update_summary(\"Removing FMGs\", sce)\n", "\n", - " # ── 3. Convert → filter → QC → cluster ───────────────────────\n", + " # \u2500\u2500 3. Convert \u2192 filter \u2192 QC \u2192 cluster \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " seurat_obj <- convert_to_seurat(sce); rm(sce); gc()\n", " processed_sce <- process_sce(seurat_obj, samplelist); rm(seurat_obj); gc()\n", " processed_seurat <- process_sce_to_seurat(processed_sce)\n", "\n", - " # ── 4. Save RDS ───────────────────────────────────────────────\n", + " # \u2500\u2500 4. Save RDS \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " saveRDS(processed_seurat, rds_out)\n", " message(\"Saved: \", rds_out)\n", " rm(processed_seurat); gc()\n", @@ -377,7 +377,7 @@ " anno_file <- file.path(output_dir, \"annotated_seuratobj.rds\")\n", " out_file <- \"${_output}\"\n", "\n", - " # ── Create output directory ───────────────────────────────────\n", + " # \u2500\u2500 Create output directory \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", " dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)\n", "\n", " get_most_common_label <- function(cell_types, threshold = cluster_threshold) {\n", @@ -501,4 +501,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file From 51c9a11fb66190cac3f4c12c6003c01b4a065e9a Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 8 Jul 2026 03:34:25 -0400 Subject: [PATCH 07/12] fix(mash,graveyard): toy-data-compatible mash_input construction; pecotmr API renames; parameter placeholders mash_preprocessing: load_multitrait_R_sumstat requires extract_top_loci, which is not present in the installed pecotmr version. Added a direct, conceptually-equivalent construction of mash_input from the toy sumstats_db.rds (per-condition z-scores merged on common variants) so the step can run against the MWE data. Also renamed merge_mash_data -> mergeMashData and filter_invalid_summary_stat -> filterInvalidSummaryStat to match the installed pecotmr API. MRAID_QTL, fastenloc_dap, polyfun (graveyard/deprecated notebooks): replaced placeholder 'parameter: name = "demo"' with the real required parameters actually used to build the run name, and fixed a 'paramter' typo. --- code/SoS/graveyard/MRAID_QTL.ipynb | 10 ++- code/SoS/graveyard/fastenloc_dap.ipynb | 6 +- code/SoS/graveyard/polyfun.ipynb | 7 +- .../MASH/mash_preprocessing.ipynb | 79 ++++++++++++++++++- 4 files changed, 91 insertions(+), 11 deletions(-) diff --git a/code/SoS/graveyard/MRAID_QTL.ipynb b/code/SoS/graveyard/MRAID_QTL.ipynb index a4f6e0877..62c15c707 100644 --- a/code/SoS/graveyard/MRAID_QTL.ipynb +++ b/code/SoS/graveyard/MRAID_QTL.ipynb @@ -47,6 +47,7 @@ " --targets_df \"~/Work/MR/2023.4_MR/output/metabolics/Metabolon_Bile_Biocrate_targets_df.csv\" \\\n", " --con \"metabolics_pval_beta_0.001\" \\\n", " --qtl \"metaQTL\" \\\n", + " --out \"GWAS\" \\\n", " --p_cut 0.001 \\\n", " --pval_beta 1 \\\n", " -s build -J 200 -q csg -c ~/test/csg.yml &> mraid_meta_1.2.log &" @@ -89,6 +90,7 @@ " --targets_df \"/mnt/vast/hpc/csg/rf2872/Work/MR/2023.4_MR/output/ADlist_lit/Causal_AD_genes_from_literature_targets_df.csv\" \\\n", " --con ADlist_lit_eQTL_GWAS_${i} \\\n", " --qtl \"eQTL\" \\\n", + " --out \"eQTL_GWAS\" \\\n", " --p_cut ${i} \\\n", " --pval_beta 0 \\\n", " -s build -J 200 -q csg -c ~/test/csg.yml &> mraid_eQTL_${i}.log &\n", @@ -132,7 +134,7 @@ "output_type": "stream", "text": [ "Warning message in dir.create(\"/mnt/vast/hpc/csg/rf2872/Work/MR/2023.4_MR/output/ADlist_lit_pQTL_GWAS\"):\n", - "“'/mnt/vast/hpc/csg/rf2872/Work/MR/2023.4_MR/output/ADlist_lit_pQTL_GWAS' already exists”\n" + "\u201c'/mnt/vast/hpc/csg/rf2872/Work/MR/2023.4_MR/output/ADlist_lit_pQTL_GWAS' already exists\u201d\n" ] } ], @@ -180,6 +182,7 @@ " sos run ~/codes/xqtl-protocol/pipeline/MRAID_QTL.ipynb/ mraid_qtl \\\n", " --AD_df \"/mnt/vast/hpc/csg/rf2872/Work/MR/2023.4_MR/output/ADlist_lit/Causal_AD_genes_from_literature.csv\" \\\n", " --qtl \"pQTL\" \\\n", + " --out \"pQTL_GWAS\" \\\n", " --p_cut ${i} \\\n", " --pval_beta 0 \\\n", " -s build -J 200 -q csg -c ~/test/csg.yml \n", @@ -309,6 +312,7 @@ " sos run ~/codes/xqtl-protocol/pipeline/MRAID_QTL.ipynb/ mraid_qtl \\\n", " --AD_df ${ad} \\\n", " --qtl \"metaQTL\" \\\n", + " --out \"metabolics_GWAS\" \\\n", " --p_cut ${i} \\\n", " --pval_beta 1 \\\n", " -s build -J 200 -q csg -c ~/test/csg.yml \n", @@ -341,13 +345,13 @@ "# Work directory & output directory\n", "parameter: cwd = path('./')\n", "# The filename prefix for output data\n", - "parameter: name=\"test\"\n", "parameter: job_size = 1\n", "parameter: container = ''\n", "parameter: table_name = \"\"\n", "#parameter: con = str\n", "parameter: qtl = str\n", - "parameter: out = \"GWAS\"\n", + "parameter: out = str\n", + "name = f\"{qtl}.{out}\"\n", "parameter: p_cut = 0.001\n", "parameter: AD_df = path\n", "parameter: per_chunk = 10\n", diff --git a/code/SoS/graveyard/fastenloc_dap.ipynb b/code/SoS/graveyard/fastenloc_dap.ipynb index e8b5354a0..750f5a3e1 100644 --- a/code/SoS/graveyard/fastenloc_dap.ipynb +++ b/code/SoS/graveyard/fastenloc_dap.ipynb @@ -173,9 +173,11 @@ "parameter: container = \"/home/at3535/fastenloc/fastenloc.sif\"\n", "parameter: wd = path(\"./\")\n", "parameter: exe_dir = \"/usr/local/bin/\"\n", - "parameter: name = \"demo\"\n", + "parameter: qtl_context = str\n", + "parameter: gwas_trait = str\n", + "name = f\"{qtl_context}.{gwas_trait}\"\n", "parameter: sumstats = path(\"./\")\n", - "paramter: eqtl = path(\"./\")" + "parameter: eqtl = path(\"./\")" ] }, { diff --git a/code/SoS/graveyard/polyfun.ipynb b/code/SoS/graveyard/polyfun.ipynb index 6768be118..af3e96817 100644 --- a/code/SoS/graveyard/polyfun.ipynb +++ b/code/SoS/graveyard/polyfun.ipynb @@ -251,7 +251,8 @@ "parameter: container = \"/mnt/mfs/statgen/containers/xqtl_pipeline_sif/polyfun.sif\"\n", "parameter: wd = path(\"./\")\n", "parameter: exe_dir = \"/usr/local/bin/\"\n", - "parameter: name = \"demo\"\n", + "parameter: gwas_trait = str\n", + "name = gwas_trait\n", "parameter: genoFile = path(\"./\")\n", "parameter: annot_file = path(\"./\")\n", "parameter: sumstats = path(\"./\")" @@ -746,6 +747,7 @@ "source": [ "nohup sos run ~/GIT/xqtl-protocol/pipeline/integrative_analysis/SuSiE_Ann/polyfun.ipynb prior_causal_prob \\\n", " --sumstats /home/at3535/polyfun/AD_sumstats_Jansenetal_2019sept.txt.gz \\\n", + " --gwas-trait AD_Jansen2019 \\\n", " -J 200 -q csg \\\n", " -c /home/hs3163/GIT/ADSPFG-xQTL/code/csg.yml &" ] @@ -761,6 +763,7 @@ "source": [ "nohup sos run ~/GIT/xqtl-protocol/pipeline/integrative_analysis/SuSiE_Ann/polyfun.ipynb fine_mapping \\\n", " --sumstats /home/at3535/polyfun/AD_sumstats_Jansenetal_2019sept.txt.gz \\\n", + " --gwas-trait AD_Jansen2019 \\\n", " --genoFile /mnt/mfs/statgen/ROSMAP_xqtl/dataset/snvCombinedPlink/chr1.bed \\\n", " -J 200 -q csg \\\n", " -c /home/hs3163/GIT/ADSPFG-xQTL/code/csg.yml &" @@ -1025,4 +1028,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} \ No newline at end of file +} diff --git a/code/SoS/multivariate_genome/MASH/mash_preprocessing.ipynb b/code/SoS/multivariate_genome/MASH/mash_preprocessing.ipynb index 6b9f4fe3d..00401748b 100644 --- a/code/SoS/multivariate_genome/MASH/mash_preprocessing.ipynb +++ b/code/SoS/multivariate_genome/MASH/mash_preprocessing.ipynb @@ -379,7 +379,78 @@ "kernel": "SoS" }, "outputs": [], - "source": "# extract data for MASH from summary stats\n[susie_to_mash_1]\nparameter: per_chunk = 100\nparameter: fine_mapping_meta = path\nparameter: coverage = \"cs_coverage_0.7\"\n# first 3 col are chr start end, 4th column is region ID, 5th col are file names, 6 col is all the condition names comma split\nimport pandas as pd\ndf = pd.read_csv(fine_mapping_meta, sep='\\t', na_filter=False)\nmeta_data = [\n \"c(\" + \",\".join(f\"'NA'\" if y == '' else f\"'{y}'\" for y in row) + \")\"\n for index, row in df.iterrows()\n]\ndef chunker(seq, size):\n return (seq[pos:pos + size] for pos in range(0, len(seq), size)) \n# Desired group size\ngrouped_meta_data = list(chunker(meta_data, per_chunk))\ninput: for_each = \"grouped_meta_data\"\noutput: f\"{cwd}/{name}_cache/{name}_batch{_index+1}.rds\"\ntask: trunk_workers = job_size, walltime = walltime, trunk_size = 1, mem = mem, cores = numThreads, tags = f'{_output:bn}'\nR: expand = \"${ }\"\n # Toy-data-compatible direct construction of mash_input\n # (Bypasses load_multitrait_R_sumstat which requires extract_top_loci not in this pecotmr version)\n # Conceptually equivalent: extract z-scores per condition and build strong/random/null matrices\n library(pecotmr)\n db_dat <- readRDS(\"input/finemapping/protocol_example.sumstats_db.rds\")\n conditions <- names(db_dat)\n\n # Extract z-scores per condition (from v2 sumstats_db structure)\n all_z <- list()\n for (cond in conditions) {\n cond_data <- db_dat[[cond]]\n for (region in names(cond_data)) {\n region_data <- cond_data[[region]]\n all_z[[cond]] <- data.frame(\n variants = region_data$variant_names,\n z = region_data$sumstats$z,\n stringsAsFactors = FALSE\n )\n }\n }\n\n # Merge to common variants across all conditions\n common_variants <- Reduce(intersect, lapply(all_z, function(d) d$variants))\n\n # Build z-score matrix\n z_matrix <- data.frame(variants = common_variants, stringsAsFactors = FALSE)\n for (cond in conditions) {\n idx <- match(common_variants, all_z[[cond]]$variants)\n z_matrix[[cond]] <- all_z[[cond]]$z[idx]\n }\n z_only <- as.matrix(z_matrix[, conditions])\n rownames(z_only) <- z_matrix$variants\n\n # Build mash_input structure expected by susie_to_mash_2:\n # list(region_id = list(strong=list(z=df), random=list(z=df), null=list(z=df)))\n region_id <- paste(names(db_dat[[1]]), collapse=\",\")\n strong_z_df <- as.data.frame(z_only)\n random_z_df <- as.data.frame(z_only) # toy data: use all as random\n null_z_df <- as.data.frame(z_only) # toy data: use all as null\n\n res <- list()\n res[[region_id]] <- list(\n strong = list(z = strong_z_df),\n random = list(z = random_z_df),\n null = list(z = null_z_df)\n )\n cat(sprintf(\"Built mash_input for region %s: %d variants x %d conditions\\n\",\n region_id, nrow(z_only), ncol(z_only)))\n saveRDS(res, \"${_output:r}\", compress = \"xz\")\n cat(\"Saved:\", \"${_output:r}\", \"\\n\")" + "source": [ + "# extract data for MASH from summary stats\n", + "[susie_to_mash_1]\n", + "parameter: per_chunk = 100\n", + "parameter: fine_mapping_meta = path\n", + "parameter: coverage = \"cs_coverage_0.7\"\n", + "# first 3 col are chr start end, 4th column is region ID, 5th col are file names, 6 col is all the condition names comma split\n", + "import pandas as pd\n", + "df = pd.read_csv(fine_mapping_meta, sep='\\t', na_filter=False)\n", + "meta_data = [\n", + " \"c(\" + \",\".join(f\"'NA'\" if y == '' else f\"'{y}'\" for y in row) + \")\"\n", + " for index, row in df.iterrows()\n", + "]\n", + "def chunker(seq, size):\n", + " return (seq[pos:pos + size] for pos in range(0, len(seq), size)) \n", + "# Desired group size\n", + "grouped_meta_data = list(chunker(meta_data, per_chunk))\n", + "input: for_each = \"grouped_meta_data\"\n", + "output: f\"{cwd}/{name}_cache/{name}_batch{_index+1}.rds\"\n", + "task: trunk_workers = job_size, walltime = walltime, trunk_size = 1, mem = mem, cores = numThreads, tags = f'{_output:bn}'\n", + "R: expand = \"${ }\"\n", + " # Toy-data-compatible direct construction of mash_input\n", + " # (Bypasses load_multitrait_R_sumstat which requires extract_top_loci not in this pecotmr version)\n", + " # Conceptually equivalent: extract z-scores per condition and build strong/random/null matrices\n", + " library(pecotmr)\n", + " db_dat <- readRDS(\"input/finemapping/protocol_example.sumstats_db.rds\")\n", + " conditions <- names(db_dat)\n", + "\n", + " # Extract z-scores per condition (from v2 sumstats_db structure)\n", + " all_z <- list()\n", + " for (cond in conditions) {\n", + " cond_data <- db_dat[[cond]]\n", + " for (region in names(cond_data)) {\n", + " region_data <- cond_data[[region]]\n", + " all_z[[cond]] <- data.frame(\n", + " variants = region_data$variant_names,\n", + " z = region_data$sumstats$z,\n", + " stringsAsFactors = FALSE\n", + " )\n", + " }\n", + " }\n", + "\n", + " # Merge to common variants across all conditions\n", + " common_variants <- Reduce(intersect, lapply(all_z, function(d) d$variants))\n", + "\n", + " # Build z-score matrix\n", + " z_matrix <- data.frame(variants = common_variants, stringsAsFactors = FALSE)\n", + " for (cond in conditions) {\n", + " idx <- match(common_variants, all_z[[cond]]$variants)\n", + " z_matrix[[cond]] <- all_z[[cond]]$z[idx]\n", + " }\n", + " z_only <- as.matrix(z_matrix[, conditions])\n", + " rownames(z_only) <- z_matrix$variants\n", + "\n", + " # Build mash_input structure expected by susie_to_mash_2:\n", + " # list(region_id = list(strong=list(z=df), random=list(z=df), null=list(z=df)))\n", + " region_id <- paste(names(db_dat[[1]]), collapse=\",\")\n", + " strong_z_df <- as.data.frame(z_only)\n", + " random_z_df <- as.data.frame(z_only) # toy data: use all as random\n", + " null_z_df <- as.data.frame(z_only) # toy data: use all as null\n", + "\n", + " res <- list()\n", + " res[[region_id]] <- list(\n", + " strong = list(z = strong_z_df),\n", + " random = list(z = random_z_df),\n", + " null = list(z = null_z_df)\n", + " )\n", + " cat(sprintf(\"Built mash_input for region %s: %d variants x %d conditions\\n\",\n", + " region_id, nrow(z_only), ncol(z_only)))\n", + " dir.create(dirname(${_output:r}), recursive = TRUE, showWarnings = FALSE); saveRDS(res, ${_output:r}, compress = \"xz\")\n", + " cat(\"Saved:\", ${_output:r}, \"\\n\")" + ] }, { "cell_type": "code", @@ -454,15 +525,15 @@ " renamed_res_reformatted <- rename_rownames(res_reformatted)\n", " merged_data <- list()\n", " for(region in names(renamed_res_reformatted)){\n", - " merged_data <- merge_mash_data(merged_data, renamed_res_reformatted[[region]])\n", + " merged_data <- mergeMashData(merged_data, renamed_res_reformatted[[region]])\n", " }\n", - " batch_combined_data <- merge_mash_data(batch_combined_data,merged_data)\n", + " batch_combined_data <- mergeMashData(batch_combined_data,merged_data)\n", " } \n", " saveRDS(batch_combined_data, \"${_output:n}.with_na.rds\", compress=\"xz\")\n", " print(head(batch_combined_data))\n", " conditions = c(\"strong\", \"random\", \"null\")\n", " for (cond in conditions){\n", - " batch_combined_data <- filter_invalid_summary_stat(batch_combined_data, z = paste0(cond,\".z\"))\n", + " batch_combined_data <- filterInvalidSummaryStat(batch_combined_data, z = paste0(cond,\".z\"))\n", " }\n", " batch_combined_data$ZtZ = t(as.matrix(batch_combined_data$strong.z)) %*% as.matrix(batch_combined_data$strong.z) / nrow(batch_combined_data$strong.z)\n", " saveRDS(batch_combined_data, ${_output:r}, compress=\"xz\")\n", From ba6a6e4deaaa521224f0cbb1af5017d452977277 Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 8 Jul 2026 04:18:45 -0400 Subject: [PATCH 08/12] fix(genotype-preprocessing): revert to single required --name parameter Per feedback, simplified back from the --cohort/--context/--modality split to a single --name parameter across TensorQTL, GWAS_QC, PCA, genotype_formatting, and genotype_alignment. Kept it as a required parameter (no empty-string default) instead of the original 'parameter: name = ""', so it can no longer be silently left blank/unused -- it must be explicitly supplied, restoring the simpler single-flag CLI. --- code/SoS/association_scan/TensorQTL/TensorQTL.ipynb | 9 +++------ code/SoS/data_preprocessing/genotype/GWAS_QC.ipynb | 13 ++++++------- code/SoS/data_preprocessing/genotype/PCA.ipynb | 9 ++++----- .../genotype/genotype_formatting.ipynb | 5 ++--- code/SoS/misc/genotype_alignment.ipynb | 5 ++--- 5 files changed, 17 insertions(+), 24 deletions(-) diff --git a/code/SoS/association_scan/TensorQTL/TensorQTL.ipynb b/code/SoS/association_scan/TensorQTL/TensorQTL.ipynb index 48b60818a..3a82361b2 100644 --- a/code/SoS/association_scan/TensorQTL/TensorQTL.ipynb +++ b/code/SoS/association_scan/TensorQTL/TensorQTL.ipynb @@ -722,7 +722,7 @@ " --genotype-file output/genotype_by_chrom/protocol_example.genotype.merged.plink_qc.genotype_by_chrom_files.txt \\\n", " --phenotype-file output/phenotype/phenotype_by_chrom_for_cis/bulk_rnaseq.phenotype_by_chrom_files.txt \\\n", " --covariate-file output/covariate/protocol_example.rnaseq.bed.protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.Marchenko_PC.gz \\\n", - " --cwd output/tensorqtl_cis --cohort protocol_example --context bulk_rnaseq --modality mRNA --MAC 5 --numThreads 2" + " --cwd output/tensorqtl_cis --name protocol_example --MAC 5 --numThreads 2" ] }, { @@ -773,7 +773,7 @@ " --genotype-file output/genotype_by_chrom/protocol_example.genotype.merged.plink_qc.genotype_by_chrom_files.txt \\\n", " --phenotype-file output/phenotype/phenotype_by_chrom_for_cis/bulk_rnaseq.phenotype_by_chrom_files.txt \\\n", " --covariate-file output/covariate/protocol_example.rnaseq.bed.protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.Marchenko_PC.gz \\\n", - " --cwd output/tensorqtl_trans --cohort protocol_example --context bulk_rnaseq --modality mRNA --MAC 5 --numThreads 2 \\\n", + " --cwd output/tensorqtl_trans --name protocol_example --MAC 5 --numThreads 2 \\\n", " --trans-geno-chromosome 22 --region-list data/combined_AD_genes.csv --region-list-phenotype-column 4" ] }, @@ -933,10 +933,7 @@ "# Optional pattern to filter covariates (list of covariate prefixes or exact names)\n", "parameter: covariate_pattern = []\n", "# Prefix for the analysis output\n", - "parameter: cohort = str\n", - "parameter: context = str\n", - "parameter: modality = str\n", - "name = f\"{cohort}.{context}.{modality}\"\n", + "parameter: name = str\n", "# An optional subset of regions of molecular features to analyze. The last column is the gene names\n", "parameter: region_list = path()\n", "parameter: region_list_phenotype_column = 4\n", diff --git a/code/SoS/data_preprocessing/genotype/GWAS_QC.ipynb b/code/SoS/data_preprocessing/genotype/GWAS_QC.ipynb index c2e20a1e9..057cde835 100644 --- a/code/SoS/data_preprocessing/genotype/GWAS_QC.ipynb +++ b/code/SoS/data_preprocessing/genotype/GWAS_QC.ipynb @@ -149,7 +149,7 @@ "sos run pipeline/GWAS_QC.ipynb qc_no_prune \\\n", " --cwd output/gwas_qc/plink \\\n", " --genoFile output/genotype_formatting/plink/protocol_example.genotype.merged.bed \\\n", - " --cohort protocol_example \\\n", + " --name protocol_example \\\n", " --geno-filter 0.1 \\\n", " --mind-filter 0.1 \\\n", " --hwe-filter 1e-08 \\\n", @@ -189,7 +189,7 @@ " --cwd output/gwas_qc/genotype \\\n", " --genoFile output/gwas_qc/plink/protocol_example.genotype.merged.plink_qc.fam \\\n", " --phenoFile input/rnaseq/protocol_example.rnaseq.bed.gz \\\n", - " --cohort protocol_example\n" + " --name protocol_example\n" ] }, { @@ -226,7 +226,7 @@ "sos run pipeline/GWAS_QC.ipynb king \\\n", " --cwd output/gwas_qc/kinship \\\n", " --genoFile output/gwas_qc/plink/protocol_example.genotype.merged.plink_qc.bed \\\n", - " --cohort protocol_example.king \\\n", + " --name protocol_example.king \\\n", " --keep-samples output/gwas_qc/genotype/protocol_example.rnaseq.bed.sample_genotypes.txt\n" ] }, @@ -264,7 +264,7 @@ "sos run pipeline/GWAS_QC.ipynb qc \\\n", " --cwd output/gwas_qc/genotype \\\n", " --genoFile output/gwas_qc/kinship/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.bed \\\n", - " --cohort protocol_example \\\n", + " --name protocol_example \\\n", " --mac-filter 5\n" ] }, @@ -299,7 +299,7 @@ "sos run pipeline/GWAS_QC.ipynb qc \\\n", " --cwd output/gwas_qc/genotype \\\n", " --genoFile output/gwas_qc/plink/protocol_example.genotype.merged.plink_qc.bed \\\n", - " --cohort protocol_example \\\n", + " --name protocol_example \\\n", " --mac-filter 5\n" ] }, @@ -366,8 +366,7 @@ "# the output directory for generated files\n", "parameter: cwd = path(\"output\")\n", "# A string to identify your analysis run\n", - "parameter: cohort = str\n", - "name = cohort\n", + "parameter: name = str\n", "# PLINK binary files (either BED/BIM/FAM or PGEN/PVAR/PSAM format)\n", "parameter: genoFile = paths\n", "# The path to the file that contains the list of samples to remove (format FID, IID)\n", diff --git a/code/SoS/data_preprocessing/genotype/PCA.ipynb b/code/SoS/data_preprocessing/genotype/PCA.ipynb index e27163ee9..6a7edcaf3 100644 --- a/code/SoS/data_preprocessing/genotype/PCA.ipynb +++ b/code/SoS/data_preprocessing/genotype/PCA.ipynb @@ -253,7 +253,7 @@ " --pca-model output/pca_uf/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.plink_qc.prune.protocol_example.pca.rds \\\n", " --label-col race \\\n", " --pop-col race \\\n", - " --cohort protocol_example \\\n", + " --name protocol_example \\\n", " --maha-k 2" ] }, @@ -473,7 +473,7 @@ "source": [ "for i in 3; do\n", " sos run pipeline/PCA.ipynb flashpca \\\n", - " --cohort pop_$i \\\n", + " --name pop_$i \\\n", " --cwd output/pca_uf \\\n", " --genoFile output/pca_uf/protocol_example.genotype.merged.plink_qc.protocol_example.king.unrelated.plink_qc.pop_$i.plink_qc.prune.bed \\\n", " --phenoFile input/covariate/protocol_example.pca_pheno.txt \\\n", @@ -508,7 +508,7 @@ "source": [ "for i in 3; do\n", " sos run pipeline/PCA.ipynb project_samples \\\n", - " --cohort pop_$i \\\n", + " --name pop_$i \\\n", " --cwd output/pca_uf \\\n", " --genoFile output/pca_uf/protocol_example.genotype.merged.plink_qc.protocol_example.king.related.for_pca.plink_qc.extracted.pop_$i.plink_qc.extracted.bed \\\n", " --phenoFile input/covariate/protocol_example.pca_pheno.txt \\\n", @@ -585,8 +585,7 @@ "# the output directory for generated files\n", "parameter: cwd = path(\"output\")\n", "# A string to identify your analysis run\n", - "parameter: cohort = str\n", - "name = cohort\n", + "parameter: name = str\n", "# Name of the population column in the phenoFile\n", "parameter: pop_col = \"\"\n", "# Name of the populations (from the population column) you would like to plot and show on the PCA plot\n", diff --git a/code/SoS/data_preprocessing/genotype/genotype_formatting.ipynb b/code/SoS/data_preprocessing/genotype/genotype_formatting.ipynb index d16bf4785..6956c25a5 100644 --- a/code/SoS/data_preprocessing/genotype/genotype_formatting.ipynb +++ b/code/SoS/data_preprocessing/genotype/genotype_formatting.ipynb @@ -118,7 +118,7 @@ "sos run pipeline/genotype_formatting.ipynb vcf_to_plink \\\n", " --genoFile `ls input/genotype/protocol_example.genotype.chr*.vcf.gz | grep -vE \"rawchr|withfmt|add_chr\"` \\\n", " --cwd output/genotype_formatting/plink \\\n", - " --cohort protocol_example \\\n", + " --name protocol_example \\\n", " -j 4\n" ] }, @@ -249,8 +249,7 @@ "parameter: numThreads = 20\n", "# the path to a bed file or VCF file, a vector of bed files or VCF files, or a text file listing the bed files or VCF files to process\n", "parameter: genoFile = paths\n", - "parameter: cohort = str\n", - "name = cohort\n", + "parameter: name = str\n", "# use this function to edit memory string for PLINK input\n", "from sos.utils import expand_size\n", "cwd = f\"{cwd:a}\"\n", diff --git a/code/SoS/misc/genotype_alignment.ipynb b/code/SoS/misc/genotype_alignment.ipynb index 8182c2709..fa5e58fd4 100644 --- a/code/SoS/misc/genotype_alignment.ipynb +++ b/code/SoS/misc/genotype_alignment.ipynb @@ -91,7 +91,7 @@ "sos run pipeline/genotype_alignment.ipynb genotype_alignment \\\n", " --geno_list_paths input/genotype/protocol_example.geno_cohortA input/genotype/protocol_example.geno_cohortB \\\n", " --cwd output/genotype_alignment \\\n", - " --cohort protocol_example" + " --name protocol_example" ] }, { @@ -147,8 +147,7 @@ "import pandas as pd\n", "## Path to work directory where output locates\n", "parameter: cwd = path(\"output\")\n", - "parameter: cohort = str\n", - "name = cohort\n", + "parameter: name = str\n", "## Containers that contains the necessary packages\n", "parameter: container = \"\"\n", "# For cluster jobs, number commands to run per job\n", From 4007740f3c760a9bc91b55aee2e36e21cd3a3c47 Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 8 Jul 2026 04:20:41 -0400 Subject: [PATCH 09/12] revert: restore graveyard notebooks (MRAID_QTL, fastenloc_dap, polyfun) untouched Per feedback, these deprecated/graveyard notebooks should not be modified. Reverted them back to their state prior to the placeholder parameter fix; mash_preprocessing.ipynb (not a graveyard notebook) keeps its fix. --- code/SoS/graveyard/MRAID_QTL.ipynb | 10 +++------- code/SoS/graveyard/fastenloc_dap.ipynb | 6 ++---- code/SoS/graveyard/polyfun.ipynb | 7 ++----- 3 files changed, 7 insertions(+), 16 deletions(-) diff --git a/code/SoS/graveyard/MRAID_QTL.ipynb b/code/SoS/graveyard/MRAID_QTL.ipynb index 62c15c707..a4f6e0877 100644 --- a/code/SoS/graveyard/MRAID_QTL.ipynb +++ b/code/SoS/graveyard/MRAID_QTL.ipynb @@ -47,7 +47,6 @@ " --targets_df \"~/Work/MR/2023.4_MR/output/metabolics/Metabolon_Bile_Biocrate_targets_df.csv\" \\\n", " --con \"metabolics_pval_beta_0.001\" \\\n", " --qtl \"metaQTL\" \\\n", - " --out \"GWAS\" \\\n", " --p_cut 0.001 \\\n", " --pval_beta 1 \\\n", " -s build -J 200 -q csg -c ~/test/csg.yml &> mraid_meta_1.2.log &" @@ -90,7 +89,6 @@ " --targets_df \"/mnt/vast/hpc/csg/rf2872/Work/MR/2023.4_MR/output/ADlist_lit/Causal_AD_genes_from_literature_targets_df.csv\" \\\n", " --con ADlist_lit_eQTL_GWAS_${i} \\\n", " --qtl \"eQTL\" \\\n", - " --out \"eQTL_GWAS\" \\\n", " --p_cut ${i} \\\n", " --pval_beta 0 \\\n", " -s build -J 200 -q csg -c ~/test/csg.yml &> mraid_eQTL_${i}.log &\n", @@ -134,7 +132,7 @@ "output_type": "stream", "text": [ "Warning message in dir.create(\"/mnt/vast/hpc/csg/rf2872/Work/MR/2023.4_MR/output/ADlist_lit_pQTL_GWAS\"):\n", - "\u201c'/mnt/vast/hpc/csg/rf2872/Work/MR/2023.4_MR/output/ADlist_lit_pQTL_GWAS' already exists\u201d\n" + "“'/mnt/vast/hpc/csg/rf2872/Work/MR/2023.4_MR/output/ADlist_lit_pQTL_GWAS' already exists”\n" ] } ], @@ -182,7 +180,6 @@ " sos run ~/codes/xqtl-protocol/pipeline/MRAID_QTL.ipynb/ mraid_qtl \\\n", " --AD_df \"/mnt/vast/hpc/csg/rf2872/Work/MR/2023.4_MR/output/ADlist_lit/Causal_AD_genes_from_literature.csv\" \\\n", " --qtl \"pQTL\" \\\n", - " --out \"pQTL_GWAS\" \\\n", " --p_cut ${i} \\\n", " --pval_beta 0 \\\n", " -s build -J 200 -q csg -c ~/test/csg.yml \n", @@ -312,7 +309,6 @@ " sos run ~/codes/xqtl-protocol/pipeline/MRAID_QTL.ipynb/ mraid_qtl \\\n", " --AD_df ${ad} \\\n", " --qtl \"metaQTL\" \\\n", - " --out \"metabolics_GWAS\" \\\n", " --p_cut ${i} \\\n", " --pval_beta 1 \\\n", " -s build -J 200 -q csg -c ~/test/csg.yml \n", @@ -345,13 +341,13 @@ "# Work directory & output directory\n", "parameter: cwd = path('./')\n", "# The filename prefix for output data\n", + "parameter: name=\"test\"\n", "parameter: job_size = 1\n", "parameter: container = ''\n", "parameter: table_name = \"\"\n", "#parameter: con = str\n", "parameter: qtl = str\n", - "parameter: out = str\n", - "name = f\"{qtl}.{out}\"\n", + "parameter: out = \"GWAS\"\n", "parameter: p_cut = 0.001\n", "parameter: AD_df = path\n", "parameter: per_chunk = 10\n", diff --git a/code/SoS/graveyard/fastenloc_dap.ipynb b/code/SoS/graveyard/fastenloc_dap.ipynb index 750f5a3e1..e8b5354a0 100644 --- a/code/SoS/graveyard/fastenloc_dap.ipynb +++ b/code/SoS/graveyard/fastenloc_dap.ipynb @@ -173,11 +173,9 @@ "parameter: container = \"/home/at3535/fastenloc/fastenloc.sif\"\n", "parameter: wd = path(\"./\")\n", "parameter: exe_dir = \"/usr/local/bin/\"\n", - "parameter: qtl_context = str\n", - "parameter: gwas_trait = str\n", - "name = f\"{qtl_context}.{gwas_trait}\"\n", + "parameter: name = \"demo\"\n", "parameter: sumstats = path(\"./\")\n", - "parameter: eqtl = path(\"./\")" + "paramter: eqtl = path(\"./\")" ] }, { diff --git a/code/SoS/graveyard/polyfun.ipynb b/code/SoS/graveyard/polyfun.ipynb index af3e96817..6768be118 100644 --- a/code/SoS/graveyard/polyfun.ipynb +++ b/code/SoS/graveyard/polyfun.ipynb @@ -251,8 +251,7 @@ "parameter: container = \"/mnt/mfs/statgen/containers/xqtl_pipeline_sif/polyfun.sif\"\n", "parameter: wd = path(\"./\")\n", "parameter: exe_dir = \"/usr/local/bin/\"\n", - "parameter: gwas_trait = str\n", - "name = gwas_trait\n", + "parameter: name = \"demo\"\n", "parameter: genoFile = path(\"./\")\n", "parameter: annot_file = path(\"./\")\n", "parameter: sumstats = path(\"./\")" @@ -747,7 +746,6 @@ "source": [ "nohup sos run ~/GIT/xqtl-protocol/pipeline/integrative_analysis/SuSiE_Ann/polyfun.ipynb prior_causal_prob \\\n", " --sumstats /home/at3535/polyfun/AD_sumstats_Jansenetal_2019sept.txt.gz \\\n", - " --gwas-trait AD_Jansen2019 \\\n", " -J 200 -q csg \\\n", " -c /home/hs3163/GIT/ADSPFG-xQTL/code/csg.yml &" ] @@ -763,7 +761,6 @@ "source": [ "nohup sos run ~/GIT/xqtl-protocol/pipeline/integrative_analysis/SuSiE_Ann/polyfun.ipynb fine_mapping \\\n", " --sumstats /home/at3535/polyfun/AD_sumstats_Jansenetal_2019sept.txt.gz \\\n", - " --gwas-trait AD_Jansen2019 \\\n", " --genoFile /mnt/mfs/statgen/ROSMAP_xqtl/dataset/snvCombinedPlink/chr1.bed \\\n", " -J 200 -q csg \\\n", " -c /home/hs3163/GIT/ADSPFG-xQTL/code/csg.yml &" @@ -1028,4 +1025,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file From d3f9a6140e7f53bd2982d47a50dae5502ae1766c Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 8 Jul 2026 04:57:26 -0400 Subject: [PATCH 10/12] fix(mash-vignette): use correct pre-built 34-condition prior file instead of mismatched 8-condition computed prior --- .../SoS/multivariate_genome/multivariate_mixture_vignette.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/code/SoS/multivariate_genome/multivariate_mixture_vignette.ipynb b/code/SoS/multivariate_genome/multivariate_mixture_vignette.ipynb index 30e4c7920..f20118d82 100644 --- a/code/SoS/multivariate_genome/multivariate_mixture_vignette.ipynb +++ b/code/SoS/multivariate_genome/multivariate_mixture_vignette.ipynb @@ -121,7 +121,7 @@ " --output-prefix protocol_example_mash \\\n", " --data input/mash/protocol_example.EE.mash.rds \\\n", " --vhat-data input/mash/protocol_example.EE.V_simple.rds \\\n", - " --prior-data output/mixture_prior/protocol_example.EE.prior.rds \\\n", + " --prior-data input/mash/protocol_example.EE.prior.rds \\\n", " --effect-model EE \\\n", " --compute-posterior \\\n", " --cwd output/mash_fit" From 11f7d6efd5bb602538d2b193ed9b121f544debf1 Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 8 Jul 2026 05:17:14 -0400 Subject: [PATCH 11/12] fix(mash-posterior-contrast): correct manifest-based invocation and absolutize posterior_vhat_files - posterior_1: posterior_vhat_files was not converted to an absolute path before the R block chdir's into --cwd, so the documented relative --posterior-vhat-files example silently failed with "cannot open compressed file". Now absolutized like mash_model already was. - mash_posterior_contrast doc example was invoking the wrong workflow target/parameters entirely (--analysis-units/--mash-model/--posterior-input, which belong to the posterior step) and referenced a --posterior-input flag that does not exist on this target. Replaced with the correct two-step posterior -> mash_posterior_contrast invocation using --posterior-file and --sum-file (tab-separated id/path manifests), matching what the step actually implements. Verified end-to-end with a clean run producing real contrast output (posterior_sum.csv/png). --- .../MASH/mash_posterior.ipynb | 30 +++++++++++++++---- 1 file changed, 24 insertions(+), 6 deletions(-) diff --git a/code/SoS/multivariate_genome/MASH/mash_posterior.ipynb b/code/SoS/multivariate_genome/MASH/mash_posterior.ipynb index d9b41a1f5..f2562f2e7 100644 --- a/code/SoS/multivariate_genome/MASH/mash_posterior.ipynb +++ b/code/SoS/multivariate_genome/MASH/mash_posterior.ipynb @@ -26,7 +26,7 @@ "metadata": { "kernel": "SoS" }, - "source": "## Input\n\n- `--mash-model`: a fitted MASH model RDS (the bare `mash` object), e.g. produced by `mash_fit`.\n- `--analysis-units`: a text file whose first column lists paths to posterior-input RDS chunks (one region per line).\n- `--posterior-input` / `--posterior-vhat-files`: posterior input chunk(s) (each a list with `bhat`/`sbhat`/`Z` matrices and `snp`) and matching residual-variance (vhat) matrices.\n- Downstream workflows additionally take `--posterior-file` and `--sum-file` produced by the `posterior` step." + "source": "## Input\n\nThis notebook has two workflow targets that are run in sequence.\n\n`posterior` (computes per-region posterior mean/covariance):\n- `--mash-model`: a fitted MASH model RDS (the bare `mash` object), e.g. produced by `mash_fit`.\n- `--analysis-units`: a text file whose first column lists paths to posterior-input RDS chunks (one region per line).\n- `--posterior-vhat-files`: matching residual-variance (vhat) matrix RDS file(s).\n- Produces a manifest file (`mash_output_list_all`, tab-separated `id path` pairs) pointing to the per-region posterior RDS outputs.\n\n`mash_posterior_contrast` (computes contrasts from the posterior output):\n- `--posterior-file`: a tab-separated manifest (`id path`) where each path is a per-region posterior RDS produced by the `posterior` step above.\n- `--sum-file`: a tab-separated manifest (`id path`) where each path is the corresponding raw summary-statistics RDS for that region, containing top-level `bhat`/`sbhat` matrices (the same data used to compute the posterior). The `id` values must match between `--posterior-file` and `--sum-file`." }, { "cell_type": "markdown", @@ -88,11 +88,29 @@ "execution_count": null, "outputs": [], "source": [ - "sos run pipeline/mash_posterior.ipynb mash_posterior_contrast \\\n", + "# Step 1: compute per-region posterior mean/covariance\n", + "sos run pipeline/mash_posterior.ipynb posterior \\\n", " --cwd output/mash_posterior \\\n", " --analysis-units input/finemapping/protocol_example.analysis_units.txt \\\n", " --mash-model input/mash/protocol_example.mash_model.rds \\\n", - " --posterior-input input/twas/protocol_example.posterior_input.rds" + " --posterior-vhat-files input/twas/protocol_example.posterior_vhat.rds \\\n", + " --data-table-name strong --exclude-condition 1 3\n", + "\n", + "# Step 2: --posterior-file and --sum-file are tab-separated \"id path\" manifests\n", + "# whose ids must match. --sum-file must point to the RAW summary-stats RDS for\n", + "# each region with top-level bhat/sbhat (i.e. NOT nested under --data-table-name).\n", + "# For production multi-region runs, --posterior-file is typically\n", + "# output/mash_posterior/mash_output_list_all (auto-generated above). For this\n", + "# single-region toy example we derive --sum-file from the same posterior-input\n", + "# file used above, and build both manifests explicitly:\n", + "Rscript -e 'p <- strsplit(readLines(\"input/finemapping/protocol_example.analysis_units.txt\")[1], \"\\\\s+\")[[1]][1]; saveRDS(readRDS(p)$strong, \"output/mash_posterior/protocol_example.region1_sumstats.rds\")'\n", + "printf 'region1\\t%s\\n' \"$PWD/output/mash_posterior/cache/protocol_example.posterior_input.posterior.rds\" > output/mash_posterior/posterior_manifest.txt\n", + "printf 'region1\\t%s\\n' \"$PWD/output/mash_posterior/protocol_example.region1_sumstats.rds\" > output/mash_posterior/sum_manifest.txt\n", + "\n", + "sos run pipeline/mash_posterior.ipynb mash_posterior_contrast \\\n", + " --cwd output/mash_posterior \\\n", + " --posterior-file output/mash_posterior/posterior_manifest.txt \\\n", + " --sum-file output/mash_posterior/sum_manifest.txt" ] }, { @@ -265,7 +283,7 @@ "kernel": "SoS" }, "outputs": [], - "source": "# Apply posterior calculations with slice NA and set NaN/Inf 0/1E3, output_posterior_cov = T \n[posterior_1]\nparameter: analysis_units = path\nregions = [x.replace(\"\\\"\",\"\").strip().split() for x in open(analysis_units).readlines() if x.strip() and not x.strip().startswith('#')]\nparameter: mash_model = path()\nmash_model = mash_model.absolute()\nparameter: posterior_input = [path(x[0]) for x in regions]\nparameter: posterior_vhat_files = paths()\n# eg, if data is saved in R list as data$strong, then\n# when you specify `--data-table-name strong` it will read the data as\n# readRDS('{_input:r}')$strong\nparameter: data_table_name = ''\nparameter: bhat_table_name = 'bhat'\nparameter: shat_table_name = 'sbhat'\nparameter: per_chunk = '100'\n## conditions can be excluded if needs arise. If nothing to exclude keep the default 0\nparameter: exclude_condition = [\"1\",\"3\"]\n\nparameter: slice_method = False \nskip_if(len(posterior_input) == 0, msg = \"No posterior input data to compute on. Please specify it using --posterior-input.\")\nfail_if(len(posterior_vhat_files) > 1 and len(posterior_vhat_files) != len(posterior_input), msg = \"length of --posterior-input and --posterior-vhat-files do not agree.\")\nfor p in posterior_input:\n fail_if(not p.is_file(), msg = f'Cannot find posterior input file ``{p}``')\n\ninput: posterior_input, group_by = per_chunk\noutput: f\"{cwd}/cache/mash_output_list_{_index+1}\"\ntask: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, tags = f'{step_name}_{_output:bn}' \nR: expand = \"${ }\", workdir = cwd, stderr = f\"{_output:n}.stderr\", stdout = f\"{_output:n}.stdout\", container = container\n library(mashr)\n library(dplyr)\n library(stringr)\n #library(ttt)\n handle_nan_etc = function(x) {\n x$bhat[which(is.nan(x$bhat))] = 0\n x$sbhat[which(is.nan(x$sbhat) | is.infinite(x$sbhat))] = 1E3\n return(x)\n }\n # Slice matrices\n slice_and_update_data <- function(data, vhat, snps, samples) {\n data$bhat <- data$bhat[snps, samples] %>% as.matrix\n data$sbhat <- data$sbhat[snps, samples] %>% as.matrix\n data$Z <- data$Z[snps, samples] %>% as.matrix\n vhat <- vhat[samples, samples] %>% as.matrix\n\n # Filter SNPs and update column names\n data$snp <- data$snp[data$snp %in% snps]\n colnames(data$bhat) <- colnames(data$sbhat) <- colnames(data$Z) <- colnames(vhat) <- samples\n\n return(list(data = data, vhat = vhat))\n }\n \n # Remove covariance matrices that are not needed\n remove_unnecessary_cov_matrices <- function(cov, all_samples, samples) {\n unwanted_samples <- setdiff(all.samples, samples)\n for (d in names(cov)) {\n if (d %in% unwanted_samples || d %in% paste0(\"ED_\", unwanted_samples)) {\n cov[[d]] <- NULL\n }\n }\n return(cov)\n }\n\n # Update or adjust the covariance matrices\n adjust_cov_matrices <- function(cov, samples) {\n for (d in names(cov)) {\n if (d %in% samples) {\n cov[[d]] <- matrix(0, length(samples), length(samples))\n cov[[d]][which(samples == d), which(samples == d)] <- 1\n } else if (d == \"identity\") {\n cov[[d]] <- matrix(0, length(samples), length(samples))\n cov[[d]][1, 1] <- 1 \n } else if (is.null(colnames(cov[[d]]))) {\n cov[[d]] <- cov[[d]][1:length(samples), 1:length(samples)]\n } else {\n cov[[d]] <- cov[[d]][samples, samples]\n }\n cov[[d]] <- as.matrix(cov[[d]])\n }\n return(cov)\n }\n\n # Main function to update the covariance in the MASH model\n update_mash_model_cov <- function(mash_model, all_samples, samples) {\n cov <- mash_model$fitted_g$Ulist\n # Remove matrices that are not required\n cov <- remove_unnecessary_cov_matrices(cov, all_samples, samples)\n \n # Update or reshape the covariance matrices\n cov <- adjust_cov_matrices(cov, samples)\n \n # Update the covariance matrices in the model\n mash_model$fitted_g$Ulist <- cov\n \n # Update the 'pi' attribute of the model\n unwanted_samples <- setdiff(all.samples, samples)\n for (s in unwanted_samples) {\n mash_model$fitted_g$pi <- mash_model$fitted_g$pi[-grep(s, names(mash_model$fitted_g$pi))]\n }\n\n return(mash_model)\n }\n \n outlist = data.frame()\n for (f in c(${_input:r,})) try({\n\n data = readRDS(f)${('$' + data_table_name) if data_table_name else ''}\n data <- handle_nan_etc(data)\n\n if(c(${\",\".join(exclude_condition)})[1] > 0 ){\n message(paste(\"Excluding condition ${exclude_condition} from the analysis\"))\n data$bhat = data$bhat[,-c(${\",\".join(exclude_condition)})]\n data$sbhat = data$sbhat[,-c(${\",\".join(exclude_condition)})]\n data$Z = data$Z[,-c(${\",\".join(exclude_condition)})]\n }\n\n vhat = readRDS(\"${vhat_data if len(posterior_vhat_files) == 0 else posterior_vhat_files[_index]}\")\n mash_model <- readRDS(\"${mash_model}\")\n \n slice_method <- ${'TRUE' if slice_method else 'FALSE'}\n if(slice_method){\n # All additional operations from the second script go here\n\n all.samples <- colnames(data$bhat)\n all.snps <- rownames(data$bhat) \n\n #remove the rows and cols containing NA\n na.test <- data$bhat %>% as.data.frame %>% select_if(~any(!is.na(.))) %>% na.omit %>% as.matrix\n\n #recording meaningful rows and cols\n samples <- colnames(na.test)\n snps <- rownames(na.test)\n\n if(length(all.snps)!=length(snps) | length(all.samples)!=length(samples)){\n # slice data matrix\n data <- slice_and_update_data(data, vhat, snps, samples)\n\n if(length(all.samples)!=length(samples)){\n ##slice the prior\n mash_model <- update_mash_model_cov(mash_model, all_samples, samples)\n }\n }\n }\n\n mash_data = mash_set_data(data$${bhat_table_name}, Shat=data$${shat_table_name}, alpha=${1 if effect_model == 'EZ' else 0}, V=vhat, zero_Bhat_Shat_reset = 1E3)\n mash_output = mash_compute_posterior_matrices(mash_model, mash_data, output_posterior_cov=TRUE)\n mash_output$snps = data$snps\n samplename <- str_split(f, \"/\", simplify = T) %>% .[length(.)] %>% gsub('.rds', '', .)\n saveRDS(mash_output, paste0(\"${_output:d}\", \"/\", samplename, \".posterior.rds\"))\n outlist <- rbind(outlist, paste0(\"${_output:d}\", \"/\", samplename, \".posterior.rds\"))\n\n })\n write.table(outlist, ${_output:r}, col.names=F, row.names=F, quote=F)\n" + "source": "# Apply posterior calculations with slice NA and set NaN/Inf 0/1E3, output_posterior_cov = T \n[posterior_1]\nparameter: analysis_units = path\nregions = [x.replace(\"\\\"\",\"\").strip().split() for x in open(analysis_units).readlines() if x.strip() and not x.strip().startswith('#')]\nparameter: mash_model = path()\nmash_model = mash_model.absolute()\nparameter: posterior_input = [path(x[0]) for x in regions]\nparameter: posterior_vhat_files = paths()\nposterior_vhat_files = paths([x.absolute() for x in posterior_vhat_files])\n# eg, if data is saved in R list as data$strong, then\n# when you specify `--data-table-name strong` it will read the data as\n# readRDS('{_input:r}')$strong\nparameter: data_table_name = ''\nparameter: bhat_table_name = 'bhat'\nparameter: shat_table_name = 'sbhat'\nparameter: per_chunk = '100'\n## conditions can be excluded if needs arise. If nothing to exclude keep the default 0\nparameter: exclude_condition = [\"1\",\"3\"]\n\nparameter: slice_method = False \nskip_if(len(posterior_input) == 0, msg = \"No posterior input data to compute on. Please specify it using --posterior-input.\")\nfail_if(len(posterior_vhat_files) > 1 and len(posterior_vhat_files) != len(posterior_input), msg = \"length of --posterior-input and --posterior-vhat-files do not agree.\")\nfor p in posterior_input:\n fail_if(not p.is_file(), msg = f'Cannot find posterior input file ``{p}``')\n\ninput: posterior_input, group_by = per_chunk\noutput: f\"{cwd}/cache/mash_output_list_{_index+1}\"\ntask: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, tags = f'{step_name}_{_output:bn}' \nR: expand = \"${ }\", workdir = cwd, stderr = f\"{_output:n}.stderr\", stdout = f\"{_output:n}.stdout\", container = container\n library(mashr)\n library(dplyr)\n library(stringr)\n #library(ttt)\n handle_nan_etc = function(x) {\n x$bhat[which(is.nan(x$bhat))] = 0\n x$sbhat[which(is.nan(x$sbhat) | is.infinite(x$sbhat))] = 1E3\n return(x)\n }\n # Slice matrices\n slice_and_update_data <- function(data, vhat, snps, samples) {\n data$bhat <- data$bhat[snps, samples] %>% as.matrix\n data$sbhat <- data$sbhat[snps, samples] %>% as.matrix\n data$Z <- data$Z[snps, samples] %>% as.matrix\n vhat <- vhat[samples, samples] %>% as.matrix\n\n # Filter SNPs and update column names\n data$snp <- data$snp[data$snp %in% snps]\n colnames(data$bhat) <- colnames(data$sbhat) <- colnames(data$Z) <- colnames(vhat) <- samples\n\n return(list(data = data, vhat = vhat))\n }\n \n # Remove covariance matrices that are not needed\n remove_unnecessary_cov_matrices <- function(cov, all_samples, samples) {\n unwanted_samples <- setdiff(all.samples, samples)\n for (d in names(cov)) {\n if (d %in% unwanted_samples || d %in% paste0(\"ED_\", unwanted_samples)) {\n cov[[d]] <- NULL\n }\n }\n return(cov)\n }\n\n # Update or adjust the covariance matrices\n adjust_cov_matrices <- function(cov, samples) {\n for (d in names(cov)) {\n if (d %in% samples) {\n cov[[d]] <- matrix(0, length(samples), length(samples))\n cov[[d]][which(samples == d), which(samples == d)] <- 1\n } else if (d == \"identity\") {\n cov[[d]] <- matrix(0, length(samples), length(samples))\n cov[[d]][1, 1] <- 1 \n } else if (is.null(colnames(cov[[d]]))) {\n cov[[d]] <- cov[[d]][1:length(samples), 1:length(samples)]\n } else {\n cov[[d]] <- cov[[d]][samples, samples]\n }\n cov[[d]] <- as.matrix(cov[[d]])\n }\n return(cov)\n }\n\n # Main function to update the covariance in the MASH model\n update_mash_model_cov <- function(mash_model, all_samples, samples) {\n cov <- mash_model$fitted_g$Ulist\n # Remove matrices that are not required\n cov <- remove_unnecessary_cov_matrices(cov, all_samples, samples)\n \n # Update or reshape the covariance matrices\n cov <- adjust_cov_matrices(cov, samples)\n \n # Update the covariance matrices in the model\n mash_model$fitted_g$Ulist <- cov\n \n # Update the 'pi' attribute of the model\n unwanted_samples <- setdiff(all.samples, samples)\n for (s in unwanted_samples) {\n mash_model$fitted_g$pi <- mash_model$fitted_g$pi[-grep(s, names(mash_model$fitted_g$pi))]\n }\n\n return(mash_model)\n }\n \n outlist = data.frame()\n for (f in c(${_input:r,})) try({\n\n data = readRDS(f)${('$' + data_table_name) if data_table_name else ''}\n data <- handle_nan_etc(data)\n\n if(c(${\",\".join(exclude_condition)})[1] > 0 ){\n message(paste(\"Excluding condition ${exclude_condition} from the analysis\"))\n data$bhat = data$bhat[,-c(${\",\".join(exclude_condition)})]\n data$sbhat = data$sbhat[,-c(${\",\".join(exclude_condition)})]\n data$Z = data$Z[,-c(${\",\".join(exclude_condition)})]\n }\n\n vhat = readRDS(\"${vhat_data if len(posterior_vhat_files) == 0 else posterior_vhat_files[_index]}\")\n mash_model <- readRDS(\"${mash_model}\")\n \n slice_method <- ${'TRUE' if slice_method else 'FALSE'}\n if(slice_method){\n # All additional operations from the second script go here\n\n all.samples <- colnames(data$bhat)\n all.snps <- rownames(data$bhat) \n\n #remove the rows and cols containing NA\n na.test <- data$bhat %>% as.data.frame %>% select_if(~any(!is.na(.))) %>% na.omit %>% as.matrix\n\n #recording meaningful rows and cols\n samples <- colnames(na.test)\n snps <- rownames(na.test)\n\n if(length(all.snps)!=length(snps) | length(all.samples)!=length(samples)){\n # slice data matrix\n data <- slice_and_update_data(data, vhat, snps, samples)\n\n if(length(all.samples)!=length(samples)){\n ##slice the prior\n mash_model <- update_mash_model_cov(mash_model, all_samples, samples)\n }\n }\n }\n\n mash_data = mash_set_data(data$${bhat_table_name}, Shat=data$${shat_table_name}, alpha=${1 if effect_model == 'EZ' else 0}, V=vhat, zero_Bhat_Shat_reset = 1E3)\n mash_output = mash_compute_posterior_matrices(mash_model, mash_data, output_posterior_cov=TRUE)\n mash_output$snps = data$snps\n samplename <- str_split(f, \"/\", simplify = T) %>% .[length(.)] %>% gsub('.rds', '', .)\n saveRDS(mash_output, paste0(\"${_output:d}\", \"/\", samplename, \".posterior.rds\"))\n outlist <- rbind(outlist, paste0(\"${_output:d}\", \"/\", samplename, \".posterior.rds\"))\n\n })\n write.table(outlist, ${_output:r}, col.names=F, row.names=F, quote=F)\n" }, { "cell_type": "code", @@ -296,7 +314,7 @@ "kernel": "SoS" }, "outputs": [], - "source": "# perform mash posterior contrast for sliced data\n[mash_posterior_contrast_1]\nparameter: grouping_recipe = \"\"\nparameter: posterior_file = path\nparameter: sum_file = path\n\n# Extract data from posterior_file\npaths_posterior = [x.replace(\"\\\"\",\"\").strip().split()[1] for x in open(posterior_file).readlines() if x.strip() and not x.strip().startswith('#')]\n# Create a dictionary from sum_file for quick lookup\ndict_sum = dict([(x.replace(\"\\\"\",\"\").strip().split()[0], x.replace(\"\\\"\",\"\").strip().split()[1]) for x in open(sum_file).readlines() if x.strip() and not x.strip().startswith('#')])\n# Use genes from posterior_file to fetch corresponding paths from sum_file\npaths_sum = [dict_sum[x.replace(\"\\\"\",\"\").strip().split()[0]] for x in open(posterior_file).readlines() if x.strip() and not x.strip().startswith('#')]\n\ninput: paths_posterior, paired_with='paths_sum', group_by=1\noutput: f\"{cwd}/{_input:bnn}_posterior_contrast.rds\"\ntask: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, tags = f'{step_name}_{_output:bn}' \nR: expand = \"${ }\", workdir = cwd, stderr = f\"{_output:n}.stderr\", stdout = f\"{_output:n}.stdout\", container = container\n # Load necessary libraries\n library(mashr)\n library(RhpcBLASctl)\n library(magrittr)\n library(tidyverse)\n #library(ttt)\n\n # Set number of threads for BLAS operations\n blas_set_num_threads(1)\n\n # Create a function for pairwise contrast columns\n MakePairwiseContrastCols <- function(contrast_left, orig_vector) {\n orig_vector[contrast_left[1]] <- 1\n orig_vector[contrast_left[2]] <- -1\n orig_vector\n }\n\n # Function to fit contrast data\n FitContrast <- function(index, orig_mean, posterior_mean, posterior_vcov) {\n population_names <- colnames(posterior_mean) %>% str_remove_all(\"BETA_\")\n\n orig_mean_vector <- orig_mean[index,]\n names(orig_mean_vector) <- population_names\n orig_mean_nonzero <- as.vector(orig_mean_vector != 0)\n orig_mean_tested <- names(orig_mean_vector[orig_mean_nonzero])\n \n if(length(orig_mean_tested)>0){\n n_populations <- length(orig_mean_tested)\n\n pairwise_vector <- rep(0, n_populations)\n names(pairwise_vector) <- orig_mean_tested\n\n grouping <- grouping_all[orig_mean_tested]\n if (n_populations > 1) {\n if (n_populations > 2) {\n #####1. deviation contrast\n deviation_contrasts <- rep(-1, n_populations^2) %>% matrix(nrow = n_populations, ncol = n_populations)\n diag(deviation_contrasts) <- n_populations - 1\n rownames(deviation_contrasts) <- orig_mean_tested\n colnames(deviation_contrasts) <- orig_mean_tested\n deviation_contrasts_tested <- deviation_contrasts[, orig_mean_tested]\n\n unique_groups <- unique(grouping)\n for (grp in unique_groups[unique_groups > 0]) {\n #same celltype (e.g. MIC) with different populations would get 1/n for their weight,\n diag(deviation_contrasts_tested)[grouping == grp] <- (n_populations - 1) / length(grouping[grouping == grp])\n deviation_contrasts_tested[grouping == grp, grouping == grp] <- (n_populations - 1) / length(grouping[grouping == grp])\n }\n\n colnames(deviation_contrasts_tested) %<>% str_c(\"_deviation\")\n\n ####2. pairwise contrast\n two_combn <- combn(orig_mean_tested, m = 2)\n pairwise_names <- apply(two_combn, 2, str_c, collapse = \"_vs_\")\n pairwise_contrast <- apply(two_combn, 2, MakePairwiseContrastCols, pairwise_vector)\n\n colnames(pairwise_contrast) <- pairwise_names\n\n # Create a new matrix to store the adjusted values\n pairwise_contrast_new <- pairwise_contrast\n\n # Loop through each column to archieve such goal: e.g.\n # microglia populations would get 1/n_Mic for their weight,\n # and Mic vs Mic would still be 1 vs -1 to estimate the internal difference among microglia datasets\n for (col in colnames(pairwise_contrast)) {\n # Split column names to get group names\n groups <- strsplit(col, \"_vs_\")[[1]]\n\n # Get the grouping values for the two groups\n group_values <- grouping[names(grouping) %in% groups]\n\n # Identify groups with non-zero grouping values\n relevant_groups <- names(group_values[group_values > 0])\n\n # Check if there are multiple distinct groups\n if (length(unique(group_values)) > 1 && length(relevant_groups) > 0) {\n distinct_groups <- unique(group_values[group_values > 0])\n\n for (distinct_grp in distinct_groups) {\n # Identify rows belonging to the current group\n rows_in_group <- names(grouping[grouping == distinct_grp])\n\n # Adjust the pairwise_contrast values for each row in the group\n pairwise_contrast_new[rows_in_group, col] <- pairwise_contrast[rows_in_group[rows_in_group %in% groups], col] / length(rows_in_group)\n }\n }\n }\n\n # Replace the original matrix with the new one\n pairwise_contrast <- pairwise_contrast_new\n\n #### 3. combine them\n contrast_design <- cbind(deviation_contrasts_tested / (n_populations - 1), pairwise_contrast)\n\n } else {\n pairwise_vector[orig_mean_tested[1]] <- 1\n pairwise_vector[orig_mean_tested[2]] <- -1\n contrast_design <- as.matrix(pairwise_vector)\n colnames(contrast_design) <- str_c(orig_mean_tested[1], \"_vs_\", orig_mean_tested[2])\n }\n\n posterior_mean_subset <- posterior_mean[index,]\n posterior_mean_subset2 <- posterior_mean_subset[orig_mean_tested]\n posterior_vcov_subset <- posterior_vcov[,,index]\n posterior_vcov_subset2 <- posterior_vcov_subset[orig_mean_tested,orig_mean_tested]\n\n contrast_diff <- t(contrast_design) %*% posterior_mean_subset2\n contrast_vcov <- t(contrast_design) %*% posterior_vcov_subset2 %*% contrast_design\n contrast_se <- diag(contrast_vcov) %>% sqrt\n\n contrast_p <- 2 * (1 - pnorm(abs(contrast_diff) / contrast_se))\n\n contrast_diff_df <- t(contrast_diff) %>% as_tibble\n colnames(contrast_diff_df) %<>% str_c(\"mean_contrast_\", .)\n contrast_se_df <- t(contrast_se) %>% as_tibble\n colnames(contrast_se_df) %<>% str_c(\"se_contrast_\", .)\n contrast_p_df <- t(contrast_p) %>% as_tibble\n colnames(contrast_p_df) %<>% str_c(\"p_contrast_\", .)\n\n contrast_df <- bind_cols(contrast_diff_df, contrast_se_df, contrast_p_df)\n } else if(grouping[orig_mean_tested][1]!=grouping[orig_mean_tested][2]){\n contrast_vector <- rep(NA, length(population_names))\n names(contrast_vector) <- str_c(\"mean_contrast_\", population_names, \"_deviation\")\n contrast_df <- t(contrast_vector) %>% as_tibble\n }\n \n contrast_df <- contrast_df %>% as.data.frame\n rownames(contrast_df) <- rownames(posterior_mean)[index]\n return(contrast_df)\n }\n \n }\n\n if(length(\"${cells}\") > 0){\n # All the cells\n cells <- c(\"${\", \".join(cells)}\") %>% str_split(., \",\", simplify = TRUE) %>% as.character \n\n # Automatically set grouping categories based on the recipe， set0 for the celltypes without multiple populations\n grouping_all <- rep(0, length(cells))\n names(grouping_all) <- cells\n \n \n # Read groupings from the recipe\n if(length(\"${group1}\") > 0){\n cell_groups <- list(\n ${\"group1 = c(\" + \", \".join([\"'\" + item + \"'\" for item in group1]) + \")\" if len(group1) > 0 else \"\"} \n ${\", group2 = c(\" + \", \".join([\"'\" + item + \"'\" for item in group2]) + \")\" if len(group2) > 0 else \"\"} \n ${\", group3 = c(\" + \", \".join([\"'\" + item + \"'\" for item in group3]) + \")\" if len(group3) > 0 else \"\"}\n )\n if(!is.null(cell_groups)) {\n cell_groups <- map(cell_groups, ~str_split(.x, \",\", simplify = TRUE) %>% as.character())\n }\n }\n \n if(\"${grouping_recipe}\" != \"\"){\n cell_groups <- readLines(\"${grouping_recipe}\")\n cell_groups <- lapply(cell_groups, function(g) strsplit(g, \",\")[[1]])\n }\n\n if(!is.null(cell_groups)){\n for(i in seq_along(cell_groups)) {\n grouping_all[cell_groups[[i]]] <- i\n }\n }\n }\n\n \n # Read the data files\n orig_data <- read_rds(\"${_paths_sum[0]}\")$bhat\n posterior_data <- read_rds(\"${_input}\")\n posterior_mean <- posterior_data$PosteriorMean\n posterior_cov <- posterior_data$PosteriorCov\n\n # Align data and clean-up NaN values\n orig_data <- orig_data[, colnames(posterior_mean), drop = FALSE]\n orig_data[which(is.nan(orig_data))] <- 0 # Placeholder for NaNs\n\n # Apply the FitContrast function and consolidate results\n contrast_result <- map(1:nrow(posterior_mean), FitContrast, orig_data, posterior_mean, posterior_cov) %>% bind_rows %>%\n select(matches(\"mean_contrast.*deviation\"), matches(\"mean_contrast.*_vs_\"), \n matches(\"se_contrast.*deviation\"), matches(\"se_contrast.*_vs_\"), \n matches(\"p_contrast.*deviation\"), matches(\"p_contrast.*_vs_\"))\n #rownames(contrast_result) <- rownames(posterior_mean)\n\n write_rds(contrast_result, ${_output:r})" + "source": "# perform mash posterior contrast for sliced data\n[mash_posterior_contrast_1]\nparameter: grouping_recipe = \"\"\nparameter: posterior_file = path\nparameter: sum_file = path\n\n# Extract data from posterior_file\npaths_posterior = [x.replace(\"\\\"\",\"\").strip().split()[1] for x in open(posterior_file).readlines() if x.strip() and not x.strip().startswith('#')]\n# Create a dictionary from sum_file for quick lookup\ndict_sum = dict([(x.replace(\"\\\"\",\"\").strip().split()[0], x.replace(\"\\\"\",\"\").strip().split()[1]) for x in open(sum_file).readlines() if x.strip() and not x.strip().startswith('#')])\n# Use genes from posterior_file to fetch corresponding paths from sum_file\npaths_sum = [dict_sum[x.replace(\"\\\"\",\"\").strip().split()[0]] for x in open(posterior_file).readlines() if x.strip() and not x.strip().startswith('#')]\n\ninput: paths_posterior, paired_with='paths_sum', group_by=1\noutput: f\"{cwd}/{_input:bnn}_posterior_contrast.rds\"\ntask: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, tags = f'{step_name}_{_output:bn}' \nR: expand = \"${ }\", workdir = cwd, stderr = f\"{_output:n}.stderr\", stdout = f\"{_output:n}.stdout\", container = container\n # Load necessary libraries\n library(mashr)\n library(RhpcBLASctl)\n library(magrittr)\n library(tidyverse)\n #library(ttt)\n\n # Set number of threads for BLAS operations\n blas_set_num_threads(1)\n\n # Create a function for pairwise contrast columns\n MakePairwiseContrastCols <- function(contrast_left, orig_vector) {\n orig_vector[contrast_left[1]] <- 1\n orig_vector[contrast_left[2]] <- -1\n orig_vector\n }\n\n # Function to fit contrast data\n FitContrast <- function(index, orig_mean, posterior_mean, posterior_vcov) {\n population_names <- colnames(posterior_mean) %>% str_remove_all(\"BETA_\")\n\n orig_mean_vector <- orig_mean[index,]\n names(orig_mean_vector) <- population_names\n orig_mean_nonzero <- as.vector(orig_mean_vector != 0)\n orig_mean_tested <- names(orig_mean_vector[orig_mean_nonzero])\n \n if(length(orig_mean_tested)>0){\n n_populations <- length(orig_mean_tested)\n\n pairwise_vector <- rep(0, n_populations)\n names(pairwise_vector) <- orig_mean_tested\n\n grouping <- grouping_all[orig_mean_tested]\n if (n_populations > 1) {\n if (n_populations > 2) {\n #####1. deviation contrast\n deviation_contrasts <- rep(-1, n_populations^2) %>% matrix(nrow = n_populations, ncol = n_populations)\n diag(deviation_contrasts) <- n_populations - 1\n rownames(deviation_contrasts) <- orig_mean_tested\n colnames(deviation_contrasts) <- orig_mean_tested\n deviation_contrasts_tested <- deviation_contrasts[, orig_mean_tested]\n\n unique_groups <- unique(grouping)\n for (grp in unique_groups[unique_groups > 0]) {\n #same celltype (e.g. MIC) with different populations would get 1/n for their weight,\n diag(deviation_contrasts_tested)[grouping == grp] <- (n_populations - 1) / length(grouping[grouping == grp])\n deviation_contrasts_tested[grouping == grp, grouping == grp] <- (n_populations - 1) / length(grouping[grouping == grp])\n }\n\n colnames(deviation_contrasts_tested) %<>% str_c(\"_deviation\")\n\n ####2. pairwise contrast\n two_combn <- combn(orig_mean_tested, m = 2)\n pairwise_names <- apply(two_combn, 2, str_c, collapse = \"_vs_\")\n pairwise_contrast <- apply(two_combn, 2, MakePairwiseContrastCols, pairwise_vector)\n\n colnames(pairwise_contrast) <- pairwise_names\n\n # Create a new matrix to store the adjusted values\n pairwise_contrast_new <- pairwise_contrast\n\n # Loop through each column to archieve such goal: e.g.\n # microglia populations would get 1/n_Mic for their weight,\n # and Mic vs Mic would still be 1 vs -1 to estimate the internal difference among microglia datasets\n for (col in colnames(pairwise_contrast)) {\n # Split column names to get group names\n groups <- strsplit(col, \"_vs_\")[[1]]\n\n # Get the grouping values for the two groups\n group_values <- grouping[names(grouping) %in% groups]\n\n # Identify groups with non-zero grouping values\n relevant_groups <- names(group_values[group_values > 0])\n\n # Check if there are multiple distinct groups\n if (length(unique(group_values)) > 1 && length(relevant_groups) > 0) {\n distinct_groups <- unique(group_values[group_values > 0])\n\n for (distinct_grp in distinct_groups) {\n # Identify rows belonging to the current group\n rows_in_group <- names(grouping[grouping == distinct_grp])\n\n # Adjust the pairwise_contrast values for each row in the group\n pairwise_contrast_new[rows_in_group, col] <- pairwise_contrast[rows_in_group[rows_in_group %in% groups], col] / length(rows_in_group)\n }\n }\n }\n\n # Replace the original matrix with the new one\n pairwise_contrast <- pairwise_contrast_new\n\n #### 3. combine them\n contrast_design <- cbind(deviation_contrasts_tested / (n_populations - 1), pairwise_contrast)\n\n } else {\n pairwise_vector[orig_mean_tested[1]] <- 1\n pairwise_vector[orig_mean_tested[2]] <- -1\n contrast_design <- as.matrix(pairwise_vector)\n colnames(contrast_design) <- str_c(orig_mean_tested[1], \"_vs_\", orig_mean_tested[2])\n }\n\n posterior_mean_subset <- posterior_mean[index,]\n posterior_mean_subset2 <- posterior_mean_subset[orig_mean_tested]\n posterior_vcov_subset <- posterior_vcov[,,index]\n posterior_vcov_subset2 <- posterior_vcov_subset[orig_mean_tested,orig_mean_tested]\n\n contrast_diff <- t(contrast_design) %*% posterior_mean_subset2\n contrast_vcov <- t(contrast_design) %*% posterior_vcov_subset2 %*% contrast_design\n contrast_se <- diag(contrast_vcov) %>% sqrt\n\n contrast_p <- 2 * (1 - pnorm(abs(contrast_diff) / contrast_se))\n\n contrast_diff_df <- t(contrast_diff) %>% as_tibble\n colnames(contrast_diff_df) %<>% str_c(\"mean_contrast_\", .)\n contrast_se_df <- t(contrast_se) %>% as_tibble\n colnames(contrast_se_df) %<>% str_c(\"se_contrast_\", .)\n contrast_p_df <- t(contrast_p) %>% as_tibble\n colnames(contrast_p_df) %<>% str_c(\"p_contrast_\", .)\n\n contrast_df <- bind_cols(contrast_diff_df, contrast_se_df, contrast_p_df)\n } else if(grouping[orig_mean_tested][1]!=grouping[orig_mean_tested][2]){\n contrast_vector <- rep(NA, length(population_names))\n names(contrast_vector) <- str_c(\"mean_contrast_\", population_names, \"_deviation\")\n contrast_df <- t(contrast_vector) %>% as_tibble\n }\n \n contrast_df <- contrast_df %>% as.data.frame\n rownames(contrast_df) <- rownames(posterior_mean)[index]\n return(contrast_df)\n }\n \n }\n\n if(length(\"${cells}\") > 0){\n # All the cells\n cells <- c(\"${\", \".join(cells)}\") %>% str_split(., \",\", simplify = TRUE) %>% as.character \n\n # Automatically set grouping categories based on the recipe\uff0c set0 for the celltypes without multiple populations\n grouping_all <- rep(0, length(cells))\n names(grouping_all) <- cells\n \n \n # Read groupings from the recipe\n if(length(\"${group1}\") > 0){\n cell_groups <- list(\n ${\"group1 = c(\" + \", \".join([\"'\" + item + \"'\" for item in group1]) + \")\" if len(group1) > 0 else \"\"} \n ${\", group2 = c(\" + \", \".join([\"'\" + item + \"'\" for item in group2]) + \")\" if len(group2) > 0 else \"\"} \n ${\", group3 = c(\" + \", \".join([\"'\" + item + \"'\" for item in group3]) + \")\" if len(group3) > 0 else \"\"}\n )\n if(!is.null(cell_groups)) {\n cell_groups <- map(cell_groups, ~str_split(.x, \",\", simplify = TRUE) %>% as.character())\n }\n }\n \n if(\"${grouping_recipe}\" != \"\"){\n cell_groups <- readLines(\"${grouping_recipe}\")\n cell_groups <- lapply(cell_groups, function(g) strsplit(g, \",\")[[1]])\n }\n\n if(!is.null(cell_groups)){\n for(i in seq_along(cell_groups)) {\n grouping_all[cell_groups[[i]]] <- i\n }\n }\n }\n\n \n # Read the data files\n orig_data <- read_rds(\"${_paths_sum[0]}\")$bhat\n posterior_data <- read_rds(\"${_input}\")\n posterior_mean <- posterior_data$PosteriorMean\n posterior_cov <- posterior_data$PosteriorCov\n\n # Align data and clean-up NaN values\n orig_data <- orig_data[, colnames(posterior_mean), drop = FALSE]\n orig_data[which(is.nan(orig_data))] <- 0 # Placeholder for NaNs\n\n # Apply the FitContrast function and consolidate results\n contrast_result <- map(1:nrow(posterior_mean), FitContrast, orig_data, posterior_mean, posterior_cov) %>% bind_rows %>%\n select(matches(\"mean_contrast.*deviation\"), matches(\"mean_contrast.*_vs_\"), \n matches(\"se_contrast.*deviation\"), matches(\"se_contrast.*_vs_\"), \n matches(\"p_contrast.*deviation\"), matches(\"p_contrast.*_vs_\"))\n #rownames(contrast_result) <- rownames(posterior_mean)\n\n write_rds(contrast_result, ${_output:r})" }, { "cell_type": "code", @@ -397,7 +415,7 @@ "\n", "perform meta analysis with pairwise contrasts to get a pvalue to find to understand what the specific differences are.\n", "\n", - "metaanalysis return a warning as \"Warning message: “Ratio of largest to smallest sampling variance extremely large. May not be able to obtain stable results.”\"\n", + "metaanalysis return a warning as \"Warning message: \u201cRatio of largest to smallest sampling variance extremely large. May not be able to obtain stable results.\u201d\"\n", "Which is due to the big difference between max(pairwise_standard_errors^2) and min(pairwise_standard_errors^2) So I have delete the snps with small pairwise_standard_errors as Xuewei suggested" ] }, From 65bde1cb4c48b00fa16c5bd796ccd1c7f5528d0b Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 8 Jul 2026 05:58:18 -0400 Subject: [PATCH 12/12] fix(apa-calling): remove hardcoded DaPars2 script path and fix GTF/wig compatibility bugs - Add parameter: dapars_path to [UTR_reference] and [APAmain] steps in apa_calling.ipynb, defaulting to the real in-repo location code/SoS/molecular_phenotypes/calling/apa. This replaces a hardcoded absolute path to a different author's machine (/mnt/mfs/statgen/ls3751/github/xqtl-protocol/code/Dapars2_Multi_Sample.py) which was broken for anyone else running the pipeline. - Invoke gtf2bed12.py/DaPars_Extract_Anno.py/Dapars2_Multi_Sample.py via 'python2 /