Imputes missing values in a molecular phenotype matrix, so downstream models see a complete matrix.
Overview¶
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.
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. 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.
When to run it. After phenotype QC and normalisation, before covariate preprocessing and the association scan.
Input¶
--phenoFile: the molecular phenotype matrix in BED layout, the first four columns being#chr,start,endandIDand the rest one column per sample, with missing entries written asNA. Exampletests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gz, a toy matrix with roughly 10% of values set to missing so the methods have gaps to fill:#chr start end ID SAMPLE_001 SAMPLE_002 chr12 752578 752579 chr12_ENSG00000060237_Q9H4A3 0.238966360190167 -0.611171227886468 chr12 990508 990509 chr12_ENSG00000082805_Q8IUD2 NA -1.86313205860919--num_factor: the number of factors for theEBMFandgEBMFmodels. It must be smaller than the number of samples; the toy set has 60, so the examples use30.--cwd: the directory outputs are written to.--save_flash: off by default. When set,gEBMFalso writes the fitted flashier factor object.--no-qc-prior-to-impute: leaves quality control on by default. Passing it skips QC, which also means the QC matrix is not available to the output.
Output¶
{cwd}/{phenoFile}.imputed.bed.gz- the imputed matrix in the same BED layout as the input, bgzipped and tabix-indexed, with a.tbibeside it. Exampleoutput/phenotype_imputation_uf/protocol_example.protein.missing.bed.imputed.bed.gz, where the gaps shown in the input preview above are filled:#chr start end ID SAMPLE_001 SAMPLE_002 SAMPLE_003 SAMPLE_004 SAMPLE_005 SAMPLE_006 SAMP chr12 752578 752579 chr12_ENSG00000060237_Q9H4A3 0.238966360190167 -0.611171227886468 1. chr12 990508 990509 chr12_ENSG00000082805_Q8IUD2 -0.762712503460405 -1.86313205860919 1.{cwd}/{phenoFile}.gEBMF_factors.rds- the saved flashier factor object from agEBMFfit, written only when--save_flash Trueis given. It is a binary R object, so it has no text preview; load it withreadRDS. Exampleoutput/phenotype_imputation_uf/protocol_example.protein.missing.bed.gEBMF_factors.rds.
Minimal Working Example¶
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.
Grouped Empirical Bayes Matrix Factorization (gEBMF)¶
Grouped Empirical Bayes Matrix Factorization. Extends EBMF by fitting factors within row groups (here, by chromosome), so structure shared across molecular features in the same group is borrowed when filling in missing values. This is the recommended default for the protocol.
Timing: TBD (on toy dataset)
sos run pipeline/phenotype_imputation.ipynb gEBMF \
--phenoFile tests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gz \
--cwd output/phenotype_imputation_uf \
--num_factor 30
Empirical Bayes Matrix Factorization (EBMF)¶
Empirical Bayes Matrix Factorization. Decomposes the phenotype matrix into a small number of latent factors with adaptive shrinkage priors, then reconstructs the matrix to impute missing entries. Captures global low-rank structure across all samples and features.
Timing: TBD (on toy dataset)
sos run pipeline/phenotype_imputation.ipynb EBMF \
--phenoFile tests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gz \
--cwd output/phenotype_imputation_uf \
--num_factor 30
missForest¶
A non-parametric, iterative random-forest imputation. Each feature with missing values is predicted in turn from the other features using a random forest, repeating until the imputed values stabilize. Handles non-linear relationships but is computationally heavier.
Timing: TBD (on toy dataset)
sos run pipeline/phenotype_imputation.ipynb missforest \
--phenoFile tests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gz \
--cwd output/phenotype_imputation_uf
sos run pipeline/phenotype_imputation.ipynb missxgboost \
--phenoFile tests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gz \
--cwd output/phenotype_imputation_uf
KNN (k-nearest-neighbours)¶
k-Nearest-Neighbours imputation. A missing entry is filled with a (distance-weighted) average of the corresponding values from the k most similar samples. Simple and fast, and works well when samples cluster into similar profiles.
Timing: TBD (on toy dataset)
sos run pipeline/phenotype_imputation.ipynb knn \
--phenoFile tests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gz \
--cwd output/phenotype_imputation_uf
SoftImpute¶
Matrix-completion via iterative soft-thresholded singular value decomposition. Assumes the phenotype matrix is approximately low-rank and recovers missing entries by repeatedly shrinking the singular values. A good linear baseline for structured data.
Timing: TBD (on toy dataset)
sos run pipeline/phenotype_imputation.ipynb soft \
--phenoFile tests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gz \
--cwd output/phenotype_imputation_uf
Mean imputation¶
Mean imputation. Each missing value is replaced by the mean of the observed values for that feature. The simplest baseline; fast but ignores correlations between samples and features.
Timing: TBD (on toy dataset)
sos run pipeline/phenotype_imputation.ipynb mean \
--phenoFile tests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gz \
--cwd output/phenotype_imputation_uf
LOD (limit-of-detection) imputation¶
Limit-of-detection imputation. Missing values are replaced with a low constant derived from the smallest observed values, appropriate when missingness is driven by measurements falling below an assay’s detection threshold (e.g., low-abundance proteins/metabolites).
Timing: TBD (on toy dataset)
sos run pipeline/phenotype_imputation.ipynb lod \
--phenoFile tests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gz \
--cwd output/phenotype_imputation_uf
Drop features that are too sparse to impute¶
bed_filter_na removes rows whose missing fraction is above the threshold, so that the imputation methods are not asked to invent values for features with almost no data. Run it before a method rather than instead of one.
Timing: TBD (on toy dataset)
sos run pipeline/phenotype_imputation.ipynb bed_filter_na \
--cwd output/phenotype_imputation_uf \
--phenoFile tests/fixtures/phenotype_imputation/protocol_example.protein.missing.bed.gzCommand Interface¶
sos run pipeline/phenotype_imputation.ipynb -husage: sos run pipeline/phenotype_imputation.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:
EBMF
gEBMF
missforest
knn
soft
mean
lod
bed_filter_na
Global Workflow Options:
--modular-script-dir code/script (as path)
--cwd output (as path)
Work directory & output directory
--phenoFile VAL (as path, required)
Molecular phenotype matrix
--[no-]qc-prior-to-impute (default to True)
QC before imputation to remove unqualified features
--job-size 1 (as int)
For cluster jobs, number commands to run per job
--walltime 72h
Wall clock time expected
--mem 16G
Memory expected
--numThreads 20 (as int)
Number of threads
--container ''
--entrypoint ''
Sections
EBMF:
Workflow Options:
--prior 'ebnm_point_laplace'
prior distribution of loadings and factors
--varType 1
--num-factor 60
gEBMF:
Workflow Options:
--nCores 1
--num-factor 60
--backfit-iter 3
--[no-]save-flash (default to True)
--[no-]null-check (default to False)
missforest:
knn:
soft:
mean:
lod:
bed_filter_na:
Workflow Options:
--rank-max 50 (as int)
--lambda-hyp 30 (as int)
--impute-method soft
--tol-missing 0.05 (as float)
Tolerance of missingness rows with missing rate larger
than tol_missing will be removed, with missing rate
smaller than tol_missing will be mean_imputed. Say if we
want to keep rows with less than 5% missing, then we use
0.05 as tol_missing.Workflow implementation¶
[global]
parameter: modular_script_dir = path('code/script') # override with --modular-script-dir
# Work directory & output directory
parameter: cwd = path("output")
# Molecular phenotype matrix
parameter: phenoFile = path
# QC before imputation to remove unqualified features
parameter: qc_prior_to_impute = True
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "72h"
# Memory expected
parameter: mem = "16G"
# Number of threads
parameter: numThreads = 20
parameter: container = ""
import re
parameter: entrypoint= ""[EBMF]
# prior distribution of loadings and factors
parameter: prior = "ebnm_point_laplace"
parameter: varType = '1'
parameter: num_factor = '60'
input: phenoFile
output: f'{cwd:a}/{_input:bn}.imputed.bed.gz'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, tags = f'{step_name}_{_output:bn}', cores = numThreads
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_imputation.R \
--step EBMF \
--cwd "${cwd}" \
--phenoFile "${_input}" \
--output "${_output}" \
--qc-prior-to-impute ${"TRUE" if qc_prior_to_impute else "FALSE"} \
--prior "${prior}" \
--varType "${varType}" \
--num-factor ${num_factor} \
--numThreads ${numThreads}
[gEBMF]
parameter: nCores = '1'
parameter: num_factor = '60'
parameter: backfit_iter = '3'
parameter: save_flash = True
parameter: null_check = False
input: phenoFile
output: f'{cwd:a}/{_input:bn}.imputed.bed.gz'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, tags = f'{step_name}_{_output:bn}', cores = numThreads
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_imputation.R \
--step gEBMF \
--cwd "${cwd}" \
--phenoFile "${_input}" \
--output "${_output}" \
--qc-prior-to-impute ${"TRUE" if qc_prior_to_impute else "FALSE"} \
--num-factor ${num_factor} \
--nCores ${nCores} \
--backfit-iter ${backfit_iter} \
--save-flash ${"TRUE" if save_flash else "FALSE"} \
--null-check ${"TRUE" if null_check else "FALSE"} \
--numThreads ${numThreads}
[missforest]
input: phenoFile
output: f'{cwd:a}/{_input:bn}.imputed.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:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_imputation.R \
--step missforest \
--cwd "${cwd}" \
--phenoFile "${_input}" \
--output "${_output}" \
--qc-prior-to-impute ${"TRUE" if qc_prior_to_impute else "FALSE"} \
--numThreads ${numThreads}
[missxgboost]
input: phenoFile
output: f'{cwd:a}/{_input:bn}.imputed.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:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_imputation.R \
--step missxgboost \
--cwd "${cwd}" \
--phenoFile "${_input}" \
--output "${_output}" \
--qc-prior-to-impute ${"TRUE" if qc_prior_to_impute else "FALSE"} \
--numThreads ${numThreads}
[knn]
input: phenoFile
output: f'{cwd:a}/{_input:bn}.imputed.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:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_imputation.R \
--step knn \
--cwd "${cwd}" \
--phenoFile "${_input}" \
--output "${_output}" \
--qc-prior-to-impute ${"TRUE" if qc_prior_to_impute else "FALSE"} \
--numThreads ${numThreads}
[soft]
input: phenoFile
output: f'{cwd:a}/{_input:bn}.imputed.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:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_imputation.R \
--step soft \
--cwd "${cwd}" \
--phenoFile "${_input}" \
--output "${_output}" \
--qc-prior-to-impute ${"TRUE" if qc_prior_to_impute else "FALSE"} \
--numThreads ${numThreads}
[mean]
input: phenoFile
output: f'{cwd:a}/{_input:bn}.imputed.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:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_imputation.R \
--step mean \
--cwd "${cwd}" \
--phenoFile "${_input}" \
--output "${_output}" \
--qc-prior-to-impute ${"TRUE" if qc_prior_to_impute else "FALSE"} \
--numThreads ${numThreads}
[lod]
input: phenoFile
output: f'{cwd:a}/{_input:bn}.imputed.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:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_imputation.R \
--step lod \
--cwd "${cwd}" \
--phenoFile "${_input}" \
--output "${_output}" \
--qc-prior-to-impute ${"TRUE" if qc_prior_to_impute else "FALSE"} \
--numThreads ${numThreads}
[bed_filter_na]
parameter: rank_max = 50 # max rank estimated in the per-chr methyl matrix
parameter: lambda_hyp = 30 # hyper par, indicating the importance of the nuclear norm
parameter: impute_method = "soft"
# Tolerance of missingness rows with missing rate larger than tol_missing will be removed,
# with missing rate smaller than tol_missing will be mean_imputed. Say if we want to keep rows with less than 5% missing, then we use 0.05 as tol_missing.
parameter: tol_missing = 0.05
input: phenoFile
output: f'{cwd:a}/{_input:nn}.filter_na.{impute_method}.bed.gz'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
Rscript ${modular_script_dir}/data_preprocessing/phenotype/phenotype_imputation.R \
--step bed_filter_na \
--cwd "${cwd}" \
--phenoFile "${_input}" \
--output "${_output}" \
--rank-max ${rank_max} \
--lambda-hyp ${lambda_hyp} \
--impute-method "${impute_method}" \
--tol-missing ${tol_missing} \
--numThreads ${numThreads}