Infers hidden factors - latent technical and biological covariates - from a molecular phenotype matrix.
Overview¶
Molecular phenotype matrices carry structure nobody recorded: batch, cell-type composition, unmeasured technical variation. Left in, that structure inflates false positives in the association scan, because a variant that happens to track a batch looks like a QTL. This step estimates those hidden factors from the data and adds them to the known covariates, so the scan can regress them out.
Every method here follows the same shape: residualise the phenotype on the covariates you already have, estimate factors from what is left, then stack the chosen factors back on top of the known covariates. They differ only in how factors are estimated and how many are kept.
When to run it. After known covariates are merged and before the association scan, on the phenotype matrix the scan will use.
Input¶
--phenoFile: the molecular phenotype matrix in UCSCbed.gzform, one row per feature and one column per sample. Exampletests/fixtures/phenotype_formatting/protocol_example.rnaseq.bed.bed.gz.--covFile: the merged known-covariate matrix, an#idheader plus one column per sample, as written bycovariate_formatting. Exampleoutput/covariate/protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.gz:#id SAMPLE_001 SAMPLE_002 SAMPLE_003 sex 1 1 1 age 90.97 80.24 83.9 PC1 -0.4525795 -0.0450112 0.0704428The covariate file is optional for
PEERbut recommended, so that a proper xQTL model is built.--N: how many hidden factors to keep.0, the default, lets the method choose:Marchenko_PCkeeps components above the Marchenko-Pastur bound andPCAuses the Buja and Eyuboglu permutation. Give a positive number to fix the count instead, whichPEERneeds.--choose_k_method:Marchenkoor the permutation alternative, selecting which rule the shared PCA step applies. It also names the output file.--mean-impute-missing: off by default. When set, missing phenotype values are mean-imputed rather than the feature being dropped.--cwd: the directory outputs are written to. Must be a full path.
Output¶
{cwd}/{name}.{choose_k_method}_PC.gz- the deliverable: the known covariates with the inferredHidden_Factor_PC*rows stacked underneath, ready for the association scan.{name}is the phenotype file stem joined to the covariate file stem. Example...prune.pca.Marchenko_PC.gzhas 60 columns and 26 rows, being the header plus 2 known covariates plus the 15 genotype PCs plus the 8 hidden factors this method kept:#id SAMPLE_001 SAMPLE_002 SAMPLE_003 sex 1 1 1 age 90.97 80.24 83.9 PC1 -0.4525795337766864 -0.04501122972985456 0.07044283190889761 PC2 0.04906817968399002 -0.1674973221595418 -0.0041635808141029 ... Hidden_Factor_PC8 ... ... ...{cwd}/{name}.residual.bed.gz- the residualised phenotype matrix produced by the shared first step and consumed by whichever method follows, written with a.tbiindex.{cwd}/{name}.residual.PEER_MODEL.hd5- the fitted PEER model, fromPEERonly.{cwd}/{name}.residual.PEER.gzand{cwd}/{name}.residual.PEER.diag.pdf- the PEER factors and the convergence diagnostic plot.
Minimal Working Example¶
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.
PCA with Marchenko-Pastur (recommended)¶
The default choice. Marchenko_PC first residualises the phenotype on the known covariates, runs PCA on the residual matrix, and keeps the components whose eigenvalues exceed the Marchenko-Pastur upper bound, which is the threshold random noise alone would produce. The kept components are stacked back onto the known covariates.
Timing: TBD (on toy dataset)
sos run pipeline/covariate_hidden_factor.ipynb Marchenko_PC \
--cwd output/covariate \
--phenoFile tests/fixtures/phenotype_formatting/protocol_example.rnaseq.bed.bed.gz \
--covFile output/covariate/protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.gz \
--mean-impute-missingPEER (GTEx-style)¶
PEER infers hidden factors using a MOFA-based probabilistic model (Stegle et al.). The number of factors follows GTEx recommendations based on sample size, or can be fixed with --N. This method runs natively here (it uses the mofapy2 Python package), so the example below is fully runnable on the toy data and produces the factor matrix plus a convergence-diagnostic PDF.
Timing: TBD (on toy dataset)
sos run pipeline/covariate_hidden_factor.ipynb PEER \
--cwd output/covariate \
--phenoFile tests/fixtures/phenotype_formatting/protocol_example.rnaseq.bed.bed.gz \
--covFile output/covariate/protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.gz \
--N 3PCA with Buja & Eyuboglu permutation¶
This is the same PCA workflow as above (PCA), but the number of components to keep is chosen by the Buja & Eyuboglu permutation procedure instead of the Marchenko-Pastur threshold. Set --choose_k_method Buja_Eyuboglu. This route uses jackstraw::permutationPA (B=100 permutations), so it runs a little longer than the Marchenko variant but is fully runnable on the toy data.
Timing: TBD (on toy dataset)
sos run pipeline/covariate_hidden_factor.ipynb PCA \
--cwd output/covariate \
--phenoFile tests/fixtures/phenotype_formatting/protocol_example.rnaseq.bed.bed.gz \
--covFile output/covariate/protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.gz \
--choose_k_method Marchenko \
--mean-impute-missingCommand Interface¶
sos run pipeline/covariate_hidden_factor.ipynb -husage: sos run pipeline/covariate_hidden_factor.ipynb
[workflow_name | -t targets] [options] [workflow_options]
workflow_name: Single or combined workflows defined in this script
targets: One or more targets to generate
options: Single-hyphen sos parameters (see "sos run -h" for details)
workflow_options: Double-hyphen workflow-specific parameters
Workflows:
Marchenko_PC
PCA
PEER
Global Workflow Options:
--modular-script-dir code/script (as path)
--cwd output (as path)
The output directory for generated files. MUST BE FULL
PATH
--covFile VAL (as path, required)
Merged Covariates File
--phenoFile VAL (as path, required)
Path to the input molecular phenotype data.
--name f'{phenoFile:bnn}.{covFile:bn}'
--job-size 1 (as int)
For cluster jobs, number commands to run per job
--walltime 5h
Wall clock time expected
--mem 16G
Memory expected
--numThreads 8 (as int)
Number of threads
--container ''
Software container option
--entrypoint ''
Sections
*_1:
Workflow Options:
--[no-]mean-impute-missing (default to False)
Marchenko_PC_2, PCA_2:
Workflow Options:
--choose-k-method Marchenko
Marchenko or Marchenko
--N 0 (as int)
PEER_2:
Workflow Options:
--N 0 (as int)
N PEER factors, If do not specify or specified as 0,
default values suggested by GTEx (based on different
sample size) will be used
--iteration 1000 (as int)
Default values from PEER software: The number of max
iteration
--tol 0.001 (as float)
Prior parameters parameter: Alpha_a = 0.001 parameter:
Alpha_b = 0.1 parameter: Eps_a = 0.1 parameter: Eps_b =
10.0 Tolarance parameters
--[no-]r2-tol (default to False)
parameter: var_tol = 0.00001 minimum variance explained
criteria to drop factors while training
--convergence-mode fast
Convergence mode: Convergence mode for MOFAr "slow",
"medium" or "fast", corresponding to 1e-5%, 1e-4% or
1e-3% deltaELBO change.
--seed 999 (as int)
Seed for the MOFA variational initialisation, making
the fit reproducible
PEER_3:Workflow implementation¶
[global]
parameter: modular_script_dir = path('code/script') # override with --modular-script-dir
# The output directory for generated files. MUST BE FULL PATH
parameter: cwd = path("output")
# Merged Covariates File
parameter: covFile = path
# Path to the input molecular phenotype data.
parameter: phenoFile = path
parameter: name = f'{phenoFile:bnn}.{covFile:bn}'
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "5h"
# Memory expected
parameter: mem = "16G"
# Number of threads
parameter: numThreads = 8
# Software container option
parameter: container = ""
parameter: entrypoint= ""[*_1(computing residual on merged covariates)]
parameter: mean_impute_missing = False
input: phenoFile, covFile
output: f'{cwd}/{name}.residual.bed.gz'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, tags = f'{step_name}_{_output:bn}'
bash: expand= "${ }", stderr = f'{_output:nn}.stderr', stdout = f'{_output:nn}.stdout', container = container, entrypoint = entrypoint
Rscript ${modular_script_dir}/data_preprocessing/covariate/covariate_hidden_factor.R \
--step compute_residual \
--cwd "${cwd}" \
--phenoFile "${phenoFile}" \
--covFile "${covFile}" \
--output "${_output}" \
${"--mean-impute-missing" if mean_impute_missing else ""} \
--numThreads ${numThreads}
Principal Components Analysis on molecular phenotype matrix¶
[Marchenko_PC_2, PCA_2]
# Marchenko or Marchenko
parameter: choose_k_method = "Marchenko"
parameter: N = 0
output: f'{cwd}/{_input:bnnn}.{choose_k_method}_PC.gz'
task: trunk_workers = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
Rscript ${modular_script_dir}/data_preprocessing/covariate/covariate_hidden_factor.R \
--step Marchenko_PC \
--cwd "${cwd}" \
--residFile "${_input}" \
--covFile "${covFile}" \
--choose-k-method "${choose_k_method}" \
--output "${_output}" \
--N ${N} \
--numThreads ${numThreads}
PEER Method¶
[PEER_2]
# N PEER factors, If do not specify or specified as 0, default values suggested by
# GTEx (based on different sample size) will be used
parameter: N = 0
# Default values from PEER software:
## The number of max iteration
parameter: iteration = 1000
### Prior parameters
#parameter: Alpha_a = 0.001
#parameter: Alpha_b = 0.1
#parameter: Eps_a = 0.1
#parameter: Eps_b = 10.0
# Tolarance parameters
parameter: tol = 0.001
#parameter: var_tol = 0.00001
# minimum variance explained criteria to drop factors while training
parameter: r2_tol = False
# Convergence mode: Convergence mode for MOFAr "slow", "medium" or "fast", corresponding to 1e-5%, 1e-4% or 1e-3% deltaELBO change.
parameter: convergence_mode = "fast"
# Seed for the MOFA variational initialisation, making the fit reproducible
parameter: seed = 999
# input is the residual file from the first step
output: f'{cwd:a}/{_input:bnn}.PEER_MODEL.hd5'
task: trunk_workers = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
bash: expand= "${ }", stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout', container = container, entrypoint = entrypoint
python3 ${modular_script_dir}/data_preprocessing/covariate/covariate_hidden_factor_peer.py \
--resid-file "${_input}" \
--model-file "${_output}" \
--factors-out "${cwd:a}/${_input:bnn}.PEER.factors.tsv" \
--weights-out "${cwd:a}/${_input:bnn}.PEER.weights.tsv" \
--variance-out "${cwd:a}/${_input:bnn}.PEER.variance.tsv" \
--num-factor ${N} \
--iteration ${iteration} \
--convergence-mode "${convergence_mode}" \
--num-threads ${numThreads} \
--tol ${tol} \
--r2-tol ${r2_tol} \
--seed ${seed}
[PEER_3]
output: f'{cwd:a}/{_input:bnn}.PEER.gz', f'{cwd:a}/{_input:bnn}.PEER.diag.pdf'
task: trunk_workers = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container = container, entrypoint = entrypoint
Rscript ${modular_script_dir}/data_preprocessing/covariate/covariate_hidden_factor.R \
--step PEER_extract \
--cwd "${cwd}" \
--modelFile "${_input}" \
--covFile "${covFile}" \
--numThreads ${numThreads}