add output test fixtures

code/SoS/data_preprocessing/covariate/covariate_hidden_factor.ipynb

@@ -45,7 +45,7 @@
     "  PC1  -0.4525795  -0.0450112  0.0704428\n",
     "  ```\n",
     "\n",
-    "  The covariate file is optional for `PEER` but recommended, so that a proper xQTL model is built. `BiCV` additionally requires an indexed VCF, which the workflow generates for itself.\n",
+    "  The covariate file is optional for `PEER` but recommended, so that a proper xQTL model is built.\n",
     "\n",
     "- `--N`: how many hidden factors to keep. `0`, the default, lets the method choose: `Marchenko_PC` keeps components above the Marchenko-Pastur bound and `PCA` uses the Buja and Eyuboglu permutation. Give a positive number to fix the count instead, which `PEER` needs.\n",
     "- `--choose_k_method`: `Marchenko` or the permutation alternative, selecting which rule the shared PCA step applies. It also names the output file.\n",
@@ -75,9 +75,7 @@
     "\n",
     "- `{cwd}/{name}.residual.bed.gz` - the residualised phenotype matrix produced by the shared first step and consumed by whichever method follows, written with a `.tbi` index.\n",
     "- `{cwd}/{name}.residual.PEER_MODEL.hd5` - the fitted PEER model, from `PEER` only.\n",
-    "- `{cwd}/{name}.residual.PEER.gz` and `{cwd}/{name}.residual.PEER.diag.pdf` - the PEER factors and the convergence diagnostic plot.\n",
-    "- `{cwd}/{name}.residual.bed.fake.vcf.gz` - the indexed VCF `BiCV` builds for APEX, an intermediate rather than a result.\n",
-    "- `{cwd}/{name}.residual.BiCV.gz` - the BiCV factors."
+    "- `{cwd}/{name}.residual.PEER.gz` and `{cwd}/{name}.residual.PEER.diag.pdf` - the PEER factors and the convergence diagnostic plot.\n"
    ]
   },
   {
@@ -88,7 +86,7 @@
    "source": [
     "## Minimal Working Example\n",
     "\n",
-    "Four alternative methods, each a chain that starts from the same residualisation step and then estimates hidden factors its own way. Pick one; they are not run in sequence."
+    "Three alternative methods, each a chain that starts from the same residualisation step and then estimates hidden factors its own way. Pick one; they are not run in sequence."
    ]
   },
   {
@@ -200,42 +198,6 @@
     "    --mean-impute-missing"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "kernel": "SoS"
-   },
-   "source": [
-    "### BiCV factor analysis with APEX\n",
-    "\n",
-    "BiCV chooses the number of hidden factors by bi-cross-validation using the APEX software (Owen et al. 2016). The pipeline computes residuals on the merged covariates, builds the inputs APEX expects (an indexed VCF placeholder and a phenotype matrix), then calls `apex factor`. The number of factors follows GTEx recommendations or can be fixed with `--N`.\n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "kernel": "SoS",
-    "tags": []
-   },
-   "source": [
-    "**Timing**: TBD (on toy dataset)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "kernel": "Bash"
-   },
-   "outputs": [],
-   "source": [
-    "sos run pipeline/covariate_hidden_factor.ipynb BiCV \\\n",
-    "    --cwd output/covariate \\\n",
-    "    --phenoFile output/rnaseq/protocol_example.rnaseq.bed.bed.gz \\\n",
-    "    --covFile output/covariate/protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.gz \\\n",
-    "    --N 3"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {
@@ -279,7 +241,6 @@
     "  Marchenko_PC\n",
     "  PCA\n",
     "  PEER\n",
-    "  BiCV\n",
     "\n",
     "Global Workflow Options:\n",
     "  --modular-script-dir code/script (as path)\n",
@@ -333,15 +294,10 @@
     "                        Convergence mode: Convergence mode for MOFAr \"slow\",\n",
     "                        \"medium\" or \"fast\", corresponding to 1e-5%, 1e-4% or\n",
     "                        1e-3% deltaELBO change.\n",
+    "      --seed 999 (as int)\n",
+    "                        Seed for the MOFA variational initialisation, making\n",
+    "                        the fit reproducible\n",
     "  PEER_3:\n",
-    "  BiCV_2:\n",
-    "  BiCV_3:\n",
-    "    Workflow Options:\n",
-    "      --N 0 (as int)\n",
-    "                        N factors, if not specify, calculated based on sample\n",
-    "                        size according to GTeX\n",
-    "      --iteration 10 (as int)\n",
-    "                        The number of iteration\n",
     "```"
    ]
   },
@@ -480,11 +436,13 @@
     "parameter: r2_tol = False\n",
     "# Convergence mode: Convergence mode for MOFAr \"slow\", \"medium\" or \"fast\", corresponding to 1e-5%, 1e-4% or 1e-3% deltaELBO change.\n",
     "parameter: convergence_mode = \"fast\"\n",
+    "# Seed for the MOFA variational initialisation, making the fit reproducible\n",
+    "parameter: seed = 999\n",
     "# input is the residual file from the first step\n",
     "output: f'{cwd:a}/{_input:bnn}.PEER_MODEL.hd5'\n",
     "task: trunk_workers = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'\n",
     "bash: expand= \"${ }\", stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout', container = container, entrypoint = entrypoint\n",
-    "    python ${modular_script_dir}/data_preprocessing/covariate/covariate_hidden_factor_peer.py \\\n",
+    "    python3 ${modular_script_dir}/data_preprocessing/covariate/covariate_hidden_factor_peer.py \\\n",
     "        --resid-file \"${_input}\" \\\n",
     "        --model-file \"${_output}\" \\\n",
     "        --factors-out \"${cwd:a}/${_input:bnn}.PEER.factors.tsv\" \\\n",
@@ -495,7 +453,8 @@
     "        --convergence-mode \"${convergence_mode}\" \\\n",
     "        --num-threads ${numThreads} \\\n",
     "        --tol ${tol} \\\n",
-    "        --r2-tol ${r2_tol}\n"
+    "        --r2-tol ${r2_tol} \\\n",
+    "        --seed ${seed}\n"
    ]
   },
   {
@@ -529,78 +488,6 @@
     "* GTEx recommandation of PEER factors is [here](https://gtexportal.org/home/documentationPage)\n",
     "* Examples by PEER is at [github](https://github.com/PMBio/peer/wiki/Tutorial)"
    ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "kernel": "SoS"
-   },
-   "source": [
-    "### BiCV Factor Analysis with APEX\n",
-    "\n",
-    "The commands below were constructed and tested against the maintainers' own reading of the program, without input from the APEX authors, so treat the parameter choices as a starting point rather than as endorsed defaults."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "kernel": "SoS",
-    "vscode": {
-     "languageId": "plaintext"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "[BiCV_2]\n",
-    "output: f'{cwd:a}/{_input:bn}.fake.vcf.gz'\n",
-    "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\n",
-    "bash: expand= \"${ }\", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', container = container, entrypoint = entrypoint\n",
-    "    Rscript ${modular_script_dir}/data_preprocessing/covariate/covariate_hidden_factor.R \\\n",
-    "        --step BiCV_2 \\\n",
-    "        --cwd \"${cwd}\" \\\n",
-    "        --residFile \"${_input}\" \\\n",
-    "        --numThreads ${numThreads}"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "kernel": "SoS",
-    "vscode": {
-     "languageId": "plaintext"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "[BiCV_3]\n",
-    "# N factors, if not specify, calculated based on sample size according to GTeX\n",
-    "parameter: N = 0\n",
-    "# The number of iteration\n",
-    "parameter: iteration = 10\n",
-    "input: output_from(1), output_from(2)\n",
-    "output: f'{cwd:a}/{_input[0]:bnn}.BiCV.gz'\n",
-    "task: trunk_workers = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'\n",
-    "bash: expand= \"${ }\", stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout', container = container, entrypoint = entrypoint\n",
-    "    n_factors_file=\"${_output[0]:nn}.n_factors.txt\"\n",
-    "    Rscript ${modular_script_dir}/data_preprocessing/covariate/covariate_hidden_factor.R \\\n",
-    "        --step BiCV_nfactors \\\n",
-    "        --residFile \"${_input[0]}\" \\\n",
-    "        --N ${N} \\\n",
-    "        --output \"$n_factors_file\"\n",
-    "    n_factors=$(cat \"$n_factors_file\")\n",
-    "    out_prefix=\"${_output[0]:nn}\"\n",
-    "    apex factor \\\n",
-    "        --out \"$out_prefix\" \\\n",
-    "        --iter ${iteration} \\\n",
-    "        --factors \"$n_factors\" \\\n",
-    "        --bed \"${_input[0]}\" \\\n",
-    "        --vcf \"${_input[1]}\" \\\n",
-    "        --threads ${numThreads} \\\n",
-    "        --cov \"${covFile}\"\n",
-    "    mv -f \"$out_prefix.cov.gz\" \"${_output[0]}\""
-   ]
   }
  ],
  "metadata": {

code/SoS/data_preprocessing/phenotype/phenotype_imputation.ipynb

@@ -22,7 +22,7 @@
     "\n",
     "Molecular phenotype matrices are rarely complete: a feature can be unmeasured in a sample for technical reasons, and the models downstream generally require a full matrix. Dropping every feature with a missing entry throws away real signal, so this step imputes the gaps instead.\n",
     "\n",
-    "Empirical Bayes matrix factorisation is the primary method here, `EBMF` for a single matrix and `gEBMF` when the data are grouped, following [Qi et al., medRxiv 2023](https://doi.org/10.1101/2023.11.29.23299181). Seven alternatives are also provided - missForest, missXGBoost, k-nearest neighbours, SoftImpute, mean imputation and limit-of-detection imputation - so a method can be matched to the assay.\n",
+    "Empirical Bayes matrix factorisation is the primary method here, `EBMF` for a single matrix and `gEBMF` when the data are grouped, following [Qi et al., medRxiv 2023](https://doi.org/10.1101/2023.11.29.23299181). Six alternatives are also provided - missForest, k-nearest neighbours, SoftImpute, mean imputation and limit-of-detection imputation - so a method can be matched to the assay.\n",
     "\n",
     "**When to run it.** After phenotype QC and normalisation, before covariate preprocessing and the association scan."
    ]
@@ -75,7 +75,7 @@
    "source": [
     "## Minimal Working Example\n",
     "\n",
-    "Eight imputation methods plus a filter that drops features too sparse to impute. Pick one method. Two practical notes for the toy data: `--num_factor` for the EBMF and gEBMF factor models must be smaller than the number of samples, which is 60 here, so the examples use `--num_factor 30`; and leave quality control enabled, that is do not pass `--no-qc-prior-to-impute`, so the QC matrix used for the output is available to every method."
+    "Seven imputation methods plus a filter that drops features too sparse to impute. Pick one method. Two practical notes for the toy data: `--num_factor` for the EBMF and gEBMF factor models must be smaller than the number of samples, which is 60 here, so the examples use `--num_factor 30`; and leave quality control enabled, that is do not pass `--no-qc-prior-to-impute`, so the QC matrix used for the output is available to every method."
    ]
   },
   {
@@ -182,27 +182,6 @@
     "    --cwd output/phenotype_imputation_uf\n"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "kernel": "SoS"
-   },
-   "source": [
-    "### missXGBoost\n",
-    "\n",
-    "Iterative imputation using gradient-boosted trees (XGBoost) as the per-feature predictor instead of random forests. Often faster and more accurate than missForest on large matrices while still modeling non-linear dependencies between features.\n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "kernel": "SoS",
-    "tags": []
-   },
-   "source": [
-    "**Timing**: TBD (on toy dataset)"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -423,7 +402,6 @@
     "  EBMF\n",
     "  gEBMF\n",
     "  missforest\n",
-    "  missxgboost\n",
     "  knn\n",
     "  soft\n",
     "  mean\n",
@@ -464,7 +442,6 @@
     "      --[no-]save-flash (default to True)\n",
     "      --[no-]null-check (default to False)\n",
     "  missforest:\n",
-    "  missxgboost:\n",
     "  knn:\n",
     "  soft:\n",
     "  mean:\n",
@@ -617,7 +594,8 @@
     "output: f'{cwd:a}/{_input:bn}.imputed.bed.gz'\n",
     "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime,  mem = mem, tags = f'{step_name}_{_output:bn}'  \n",
     "bash: expand= \"${ }\", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint\n",
-    "    bash ${modular_script_dir}/data_preprocessing/phenotype/phenotype_imputation.sh missxgboost \\\n",
+    "    Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_imputation.R \\\n",
+    "        --step missxgboost \\\n",
     "        --cwd \"${cwd}\" \\\n",
     "        --phenoFile \"${_input}\" \\\n",
     "        --output \"${_output}\" \\\n",

code/SoS/misc/data_preprocessing/2_genotype_preprocessing.ipynb

@@ -18,7 +18,7 @@
     "3. **PLINK quality control** (`GWAS_QC.ipynb`): filter variants/samples by MAF, missingness (geno/mind), and Hardy-Weinberg equilibrium.\n",
     "4. **Split by chromosome** (`genotype_formatting.ipynb`): produce one PLINK file set per chromosome for downstream cis analysis.\n",
     "\n",
-    "> Note: principal-component analysis (PCA) and the genetic relationship matrix (GRM) are covered in the next mini-protocol, `3_genotype_pca.ipynb`, not here. The toy inputs use 60 samples named `SAMPLE_001 ... SAMPLE_060`, consistent with the phenotype MWE."
+    "> Note: principal-component analysis (PCA) is covered in the next mini-protocol, `3_genotype_pca.ipynb`, not here. The toy inputs use 60 samples named `SAMPLE_001 ... SAMPLE_060`, consistent with the phenotype MWE."
    ]
   },
   {

code/SoS/mnm_analysis/mnm_methods/colocboost.ipynb

@@ -517,6 +517,8 @@
     "parameter: xqtl_coloc = True\n",
     "parameter: joint_gwas = False\n",
     "parameter: separate_gwas = True\n",
+    "# RNG seed forwarded to colocboost.R (stochastic gradient boosting) for reproducible fits.\n",
+    "parameter: seed = 999\n",
     "# Summary-stats LD-mismatch QC method (none | slalom | dentist).\n",
     "parameter: qc_method = \"none\"\n",
     "# Analysis environment settings.\n",
@@ -654,6 +656,7 @@
     "        ${'' if xqtl_coloc else '--no-xqtl-coloc'} \\\n",
     "        ${'--joint-gwas' if joint_gwas else ''} \\\n",
     "        ${'--separate-gwas' if separate_gwas else ''} \\\n",
+    "        --seed ${seed} \\\n",
     "        --pip-cutoff-to-skip \"${\",\".join([str(x) for x in skip_analysis_pip_cutoff])}\" \\\n",
     "        --output ${_output}\n"
    ]