Covariate Data Formatting#
Merges known covariates, genotype principal components and hidden factors into the single file the association scan expects.
Overview#
The association scan takes a single covariate file, but covariates arrive from several
places: batch and demographic variables you supply, genotype principal components from
PCA, and hidden factors estimated by covariate_hidden_factor. This step merges them,
reconciles sample identifiers across the sources, and writes the combined matrix in the
orientation the scan expects.
When to run it. After PCA and hidden-factor estimation, immediately before QTL association testing.
Input#
--covFile: the known covariate table, samples in columns and covariates in rows, with an#idheader line. Exampleinput/covariate/protocol_example.covariates.base.tsv, showing the first 4 of 59 samples:#id SAMPLE_001 SAMPLE_002 SAMPLE_003 SAMPLE_004 sex 1 1 1 1 age 90.97 80.24 83.9 74.1
--pcaFile: the genotype PCA result written by thePCAmodule, an RDS holding the per-sample scores inpc_scores. Exampleoutput/genotype/genotype_pca/protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.rds:List of 6 $ pca_model :List of 7 ..- attr(*, "class")= chr "flashpca" $ pc_scores :'data.frame': 60 obs. of 28 variables: $ meta : chr "protocol_example.genotype.merged.plink_qc.plink_qc.prune " $ pc_cov : num [1:5, 1:5] 1.91e-02 -4.83e-18 6.02e-18 -1.35e-18 2.33e-18 ... $ pc_mean : num [1:5] 5.17e-17 -1.22e-17 -4.58e-17 7.87e-17 -3.73e-17 $ pc_median : num [1:5] 0.00271 0.00131 0.00781 -0.00328 0.00532
--k: how many genotype PCs to keep, 20 by default. The example instead reads the count off the PCA scree file, taking the last PC whose cumulative variance explained stays under 80%. Scree fileoutput/genotype/genotype_pca/protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.scree.txt:PCs PVE PVE_cum 1 0.05 0.05 2 0.05 0.1
--tol-cov: the largest fraction of missing values a covariate row may carry and still be kept, its missing entries mean-imputed rather than the row dropped. The default-1.0drops any row that has a missing value;0.4admits rows up to 40% missing.--name: the stem of the output file.--cwd: the directory the output is written to.
Output#
{cwd}/{name}.gz: the combined covariate matrix, samples in columns and covariates in rows under an#idheader, in the orientation the association scan expects. Exampleoutput/covariate/protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.gzhas 60 columns and 18 rows, the header plus the 2 supplied covariates plus the 15 PCs--kselected. First 4 columns and first 6 rows:#id SAMPLE_001 SAMPLE_002 SAMPLE_003 sex 1 1 1 age 90.97 80.24 83.9 PC1 -0.45257953377688636 -0.04501122972985456 0.07044283190889761 PC2 0.049068179683990002 -0.1674973221595418 -0.004163580814102901 PC3 -0.021287959073347167 0.011427192975949114 0.2219354934004062
Minimal Working Example#
Merge covariates and genotype PCs#
This step keeps the samples present in both the PCA result and the covariate table, then stacks the chosen genotype PCs on top of the known covariates to give the association scan a single matrix.
Timing: TBD (on toy dataset)
sos run pipeline/covariate_formatting.ipynb merge_genotype_pc \
--cwd output/covariate/ \
--pcaFile output/genotype/genotype_pca/protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.rds \
--covFile input/covariate/protocol_example.covariates.base.tsv \
--name protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca \
--tol-cov 0.4 \
--k `awk '$3 < 0.8' output/genotype/genotype_pca/protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.scree.txt | tail -1 | cut -f 1`
Command Interface#
sos run pipeline/covariate_formatting.ipynb -h
usage: sos run pipeline/covariate_formatting.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:
merge_genotype_pc
Global Workflow Options:
--modular-script-dir code/script (as path)
--cwd output (as path)
The output directory for generated files.
--covFile VAL (as path, required)
The covariate file
--job-size 1 (as int)
For cluster jobs, number commands to run per job
--walltime 5h
Wall clock time expected
--mem 2G
Memory expected
--numThreads 8 (as int)
Number of threads
--container ''
Software container option
--entrypoint ''
Sections
merge_genotype_pc:
Workflow Options:
--pcaFile VAL (as path, required)
An RDS file as the output of the genotype PCA module
--k 20 (as int)
The number of PCs to retain, by default is 20, in
practice can be the number of PC that captured more than
70% PVE
--name f'{covFile:bn}.{pcaFile:bn}'
--outliersFile . (as path)
Outliers
--[no-]remove-outliers (default to False)
--tol-cov -1.0 (as float)
Tolerance of missingness in covariates, -1 means do
nothing, otherwise for samples with covariates missing
rate larger than tol_cov will be removed, with missing
rate smaller than tol_cov will be kept.
--[no-]mean-impute (default to True)
Workflow implementation#
[global]
parameter: modular_script_dir = path('code/script') # override with --modular-script-dir
# The output directory for generated files.
parameter: cwd = path("output")
# The covariate file
parameter: covFile = path
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "5h"
# Memory expected
parameter: mem = "2G"
# Number of threads
parameter: numThreads = 8
# Software container option
parameter: container = ""
parameter: entrypoint=""
cwd = path(f"{cwd:a}")
[merge_genotype_pc]
# An RDS file as the output of the genotype PCA module
parameter: pcaFile = path
# The number of PCs to retain, by default is 20, in practice can be the number of PC that captured more than 70% PVE
parameter: k = 20
parameter: name = f'{covFile:bn}.{pcaFile:bn}'
# Outliers
parameter: outliersFile = path(".")
parameter: remove_outliers = False
# Tolerance of missingness in covariates, -1 means do nothing, otherwise for samples with covariates missing rate larger than tol_cov will be removed,
# with missing rate smaller than tol_cov will be kept.
parameter: tol_cov = -1.0
parameter: mean_impute = True
stop_if(remove_outliers and not outliersFile.is_file(), msg = "No outliers file specified, please add outliers file or remove the remove-outliers flag")
input: pcaFile, covFile
output: f'{cwd:a}/{name}.gz'
task: trunk_workers = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]: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_formatting.R \
--step merge_genotype_pc \
--cwd "${cwd}" \
--pcaFile "${pcaFile}" \
--covFile "${covFile}" \
--name "${name}" \
--k ${k} \
--outliersFile "${outliersFile}" \
${"--remove-outliers" if remove_outliers else ""} \
--tol-cov ${tol_cov} \
${"--mean-impute" if mean_impute else ""} \
--numThreads ${numThreads}