Covariate Data Formatting - FunGen-xQTL Consortium

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 #id header line. Example tests/fixtures/covariate_formatting/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 the PCA module, an RDS holding the per-sample scores in pc_scores. Example output/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 file output/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.0 drops any row that has a missing value; 0.4 admits 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 #id header, in the orientation the association scan expects. Example output/covariate/protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.gz has 60 columns and 18 rows, the header plus the 2 supplied covariates plus the 15 PCs --k selected. 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 tests/fixtures/covariate_formatting/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}