TWAS and COLOC integration (INTACT)

Run INTACT on the PTWAS and fastenloc output to integrate transcriptome-wide association (TWAS) evidence with colocalization (COLOC) evidence into a single gene-level posterior probability.

Overview

The posterior probabilities in the INTACT function are calculated through the following steps:

• Compute Bayes Factors: The function converts TWAS z-scores into Bayes factors by computing a grid of Bayes factors using different values of K (a vector of values over which Bayesian model averaging is performed, default c(1,2,4,8,16)), then performing Bayesian model averaging using the log sum exp trick.

• Compute Prior Probabilities: The prior_fun function converts GLCP into prior probabilities. prior_fun is set as linear by default.

• Compute Posterior Probabilities.

When to run it. After PTWAS and fastenloc have both been run on the same gene set, to combine their evidence.

Input

• --ptwas-file: PTWAS output, a TSV file with (at minimum) columns GENE, STAT (TWAS z-score) and SUBCLASS.

• --fastenloc-file: fastenloc output, a whitespace-delimited table with columns Gene and GLCP (gene-level colocalization probability).

• --tissue: tissue / dataset label used to name the output (e.g. DLPFC).

• --alpha: FDR significance threshold (default 0.05).

The genes are matched between the two inputs by joining on the gene name.

• --cwd: the directory outputs are written to.

Output

• {cwd}/{tissue}.INTACT.rds - the merged PTWAS and fastenloc table with the INTACT posterior probability added, sorted by descending posterior. Example output/intact/DLPFC.INTACT.rds, a data frame of 40 genes and 5 columns:

  'data.frame':  40 obs. of  5 variables:
   $ gene      : chr  "GENE021" "GENE023" "GENE004" "GENE034" ...
   $ zscore    : num  -6.55 -4.45 4.03 -3.41 -3.61 ...
   $ GLCP      : num  0.9975 0.0909 0.2677 0.8608 0.3994 ...
   $ intact_pip: num  1 0.983 0.975 0.963 0.941 ...
   $ fdr_sig   : logi  TRUE TRUE TRUE TRUE TRUE FALSE ...
  

  zscore comes from the PTWAS input and GLCP from fastenloc; intact_pip is the posterior this module computes from the two, and fdr_sig flags the genes passing --alpha.

Minimal Working Example

One workflow, intact, run once per tissue. It requires the R package INTACT from Bioconductor, provided by the pipeline environment as bioconductor-intact.

Timing: ~1-2 min (on toy dataset)

sos run pipeline/intact.ipynb intact \
    --fastenloc-file input/twas/protocol_example.fastenloc.gene.out \
    --ptwas-file input/twas/protocol_example.ptwas.output \
    --tissue DLPFC \
    --cwd output/intact

Environment note. This step requires the R package INTACT (Bioconductor), provided by the pipeline environment as bioconductor-intact. Supply real PTWAS and fastenloc output files for a production run.

Command Interface

sos run pipeline/intact.ipynb -h

usage: sos run code/SoS/pecotmr_integration/intact.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:
  intact

Global Workflow Options:
  --cwd output (as path)
                        Workdir
  --modular-script-dir code/script (as path)
                        Directory of the modular analysis scripts (code/script)
  --fastenloc-file ''
                        fastenloc output
  --ptwas-file ''
                        ptwas output
  --tissue ''
                        dataset
  --alpha 0.05 (as float)
  --QTL eQTL
                        QTL data type
  --container ''
  --job-size 1 (as int)
  --walltime 5h
  --mem 8G
  --numThreads 1 (as int)

Sections
  intact:

Workflow implementation

[global]
# Workdir
parameter: cwd = path("output")
# Directory of the modular analysis scripts (code/script)
parameter: modular_script_dir = path('code/script')
# fastenloc output
parameter: fastenloc_file = ""
# ptwas output
parameter: ptwas_file = ""
# dataset 
parameter: tissue = ''
parameter: alpha = 0.05
# QTL data type
parameter: QTL = 'eQTL'
parameter: container = ''
parameter: job_size = 1
parameter: walltime = "5h"
parameter: mem = "8G"
parameter: numThreads = 1

[intact]
input: ptwas_file, fastenloc_file
output: f'{cwd}/{tissue}.INTACT.rds'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand = "${ }", stderr = f'{_output:nn}.stderr', stdout = f'{_output:nn}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/intact.R \
        --ptwas-file ${_input[0]} \
        --fastenloc-file ${_input[1]} \
        --alpha ${alpha} \
        --output ${_output}