Hierarchical Multiple Testing

Hierarchical Multiple Testing#

Applies hierarchical multiple-testing correction to the per-gene cis-QTL results from TensorQTL and assembles the survivors into one consolidated object.

Overview#

A cis scan reports a p-value for every variant against every molecular phenotype, which is not yet a result: the variants within a gene’s window are correlated and the genes are many, so raw p-values overstate significance twice over.

The correction is hierarchical, in three steps: local adjustment of the p-values of all cis variants within each gene, global adjustment of the minimum adjusted p-value across genes, then selection of the xQTLs whose locally adjusted p-value falls below the threshold (the eigenMT-BH procedure of Huang et al. 2018, NAR 46(22):e133). The survivors are packaged as a QtlSumStats object with a regional FDR table, and the intermediate TensorQTL files are reorganised into an archive folder for book-keeping or deletion.

The same default workflow handles three flavours of upstream result: the genetic (cis) effect shown below, the interaction effect - run on *.cis_qtl_top_assoc.txt.gz with --maf-cutoff 0 --cis-window 0 and the interaction p/q-value patterns - and quantile QTL.

When to run it. Immediately after TensorQTL, before fine-mapping. Nothing downstream reads the raw scan directory directly.

Input#

The workflow is pointed at a directory rather than at files: it globs --cwd for the three TensorQTL products below using --regional-pattern, --qtl-pattern and --n-variants-suffix. output/tensorqtl_cis is one directory in that shape.

  • --cwd: the TensorQTL working directory holding the scan output. Three kinds of file are picked up from it.

    The per-gene regional result, matched by --regional-pattern (default cis_qtl[.]regional[.]tsv[.]gz$). Example output/tensorqtl_cis/protocol_example.chr22_chr22.cis_qtl.regional.tsv.gz, 21 columns and 201 rows:

    chrom  pos       n_variants  beta_shape1  beta_shape2  true_df    p_true_df             variant_id
22     11797697  6           0.9833031    5.6122484    30.858135  0.0676042590649028    chr22:11797697_C_T
22     15838003  46          0.98497593   26.038029    30.771364  0.007232871450238482  chr22:15838003_CTTTC...

    The full variant-level pairs file, matched by --qtl-pattern (default cis_qtl[.]pairs[.]tsv[.]gz$). Example output/tensorqtl_cis/protocol_example.chr22_chr22.cis_qtl.pairs.tsv.gz, 15 columns and 18319 rows:

    chrom  pos       molecular_trait_id  variant_id          tss_distance  tes_distance  af           pvalue
22     10685239  ENSG00000283047     chr22:10685239_C_T  -254148       -276098       0.059322033  0.572662489372593
22     10761227  ENSG00000283047     chr22:10761227_G_A  -178160       -200110       0.15254237   0.752156227936135

    The recounted variant tallies, matched by --n-variants-suffix (default cis_n_variants_stats[.]tsv[.]gz$). Example output/tensorqtl_cis/protocol_example.chr22_chr22.maf_0.01_window_1000000_cis_n_variants_stats.tsv.gz, 4 columns and 201 rows:

    chrom  molecular_trait_object_id  n_variants  n_variants_filtered
chr22  ENSG00000063515            44          44
chr22  ENSG00000075234            111         111

  • --output-dir: where the consolidated object and the regional FDR table are written. Required.

  • --modular-script-dir: the code/script directory holding the pecotmr_integration wrappers. Required.

  • --maf-cutoff and --cis-window: 0.01 and 1000000 by default. When either is non-zero the variant counts are recomputed from the full pairs file before Bonferroni correction, which matters because a correction taken over too wide a window, or over many low-frequency variants, overstates the number of independent tests. Set both to 0 to skip the recount and use the counts as given.

  • --pvalue-cutoff: 0.05 by default. Pairs above it are dropped before assembly to keep the object small; the per-gene minimum variant is always retained so the Bonferroni minimum stays exact.

  • --fdr-threshold: 0.05 by default, the threshold for the global significance call.

  • --study, --context, --genome: labels stamped onto the assembled QtlSumStats object, study / context / hg38 by default.

  • --af-col, --pvalue-col, --molecular-id-col: the column names to read, af / pvalue / molecular_trait_object_id by default.

  • --[no-]enable-archive: off by default. When on, the bulky TensorQTL intermediates are moved out of the working directory for book-keeping or deletion.

Output#

  • {output-dir}/{cwd basename}.qtl_association_postprocessing.rds - the consolidated QtlSumStats object holding the surviving xQTLs together with the --study, --context and --genome labels. This is the only file named in the workflow’s output: statement.

  • {output-dir}/{cwd basename}.qtl_association_postprocessing.cis_regional.fdr.tsv.gz - the per-gene regional FDR table.

Both are named from the basename of --cwd, so the example below writes tensorqtl_cis.qtl_association_postprocessing.rds and tensorqtl_cis.qtl_association_postprocessing.cis_regional.fdr.tsv.gz into output/qtlpp.

No example output ships for this step. The *.cis_regional.*.fdr.gz and *_regional_significance.* files sitting alongside the inputs in output/tensorqtl_cis were written by an earlier version of this interface and carry names the workflow no longer produces.

Minimal Working Example#

The workflow globs --cwd for its inputs rather than taking file arguments, so it is pointed straight at the TensorQTL scan directory.

Timing: TBD (on toy dataset)

sos run pipeline/qtl_association_postprocessing.ipynb default \
    --cwd output/tensorqtl_cis \
    --output-dir output/qtlpp \
    --modular-script-dir code/script \
    --maf-cutoff 0.01 --cis-window 1000000 --pvalue-cutoff 0.05 \
    --study protocol_example --context bulk_rnaseq --genome hg38

Command Interface#

sos run pipeline/qtl_association_postprocessing.ipynb -h

usage: sos run pipeline/qtl_association_postprocessing.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:
  default

Global Workflow Options:
  --cwd . (as path)
  --modular-script-dir VAL (as path, required)
                        code/script dir that holds the pecotmr_integration
                        wrappers
  --output-dir VAL (as path, required)
  --sub-dir . (as path)
  --study study
                        labels for the assembled QtlSumStats
  --context context
  --genome hg38
  --maf-cutoff 0.01 (as float)
  --cis-window 1000000 (as int)
  --af-col af
  --pvalue-col pvalue
  --molecular-id-col 'molecular_trait_object_id'
  --pvalue-cutoff 0.05 (as float)
                        drop pairs with p above this before assembling (keeps
                        the object small; the per-gene min variant is always
                        retained so the Bonferroni min is exact)
  --fdr-threshold 0.05 (as float)
  --regional-pattern 'cis_qtl[.]regional[.]tsv[.]gz$'
                        list.files() regex patterns to locate the inputs in the
                        work dir
  --qtl-pattern 'cis_qtl[.]pairs[.]tsv[.]gz$'
  --n-variants-suffix 'cis_n_variants_stats[.]tsv[.]gz$'
  --[no-]enable-archive (default to False)

Sections
  default:

Workflow implementation#

The cells below are the unmodified SoS workflow definition that the example above calls.

[global]
parameter: cwd = path(".")
# code/script dir that holds the pecotmr_integration wrappers
parameter: modular_script_dir = path
parameter: output_dir = path
parameter: sub_dir = path(".")
# labels for the assembled QtlSumStats
parameter: study = "study"
parameter: context = "context"
parameter: genome = "hg38"
parameter: maf_cutoff = 0.01
parameter: cis_window = 1000000
parameter: af_col = "af"
parameter: pvalue_col = "pvalue"
parameter: molecular_id_col = "molecular_trait_object_id"
# drop pairs with p above this before assembling (keeps the object small; the
# per-gene min variant is always retained so the Bonferroni min is exact)
parameter: pvalue_cutoff = 0.05
parameter: fdr_threshold = 0.05
# list.files() regex patterns to locate the inputs in the work dir
parameter: regional_pattern = "cis_qtl[.]regional[.]tsv[.]gz$"
parameter: qtl_pattern = "cis_qtl[.]pairs[.]tsv[.]gz$"
parameter: n_variants_suffix = "cis_n_variants_stats[.]tsv[.]gz$"
parameter: enable_archive = False

work_dir = cwd if str(sub_dir) == "." else path(f"{cwd:a}/{sub_dir}")

[default]
# Hierarchical multiple-testing correction of the cis-QTL association results,
# via pecotmr::qtlAssociationPostprocess (the wrapper reads the per-gene regional
# + per-variant pairs, assembles a per-gene QtlSumStats, and writes the enriched
# consolidated RDS + per-method regional / significant-QTL / event / summary
# tables). Replaces the old source(pecotmr/inst/code/tensorqtl_postprocessor.R).
output: f"{output_dir:a}/{cwd:b}.qtl_association_postprocessing.rds"
bash: expand = "${ }"
    Rscript ${modular_script_dir:a}/pecotmr_integration/qtl_association_postprocessing.R \
        --cwd ${work_dir:a} \
        --regional-pattern "${regional_pattern}" \
        --qtl-pattern "${qtl_pattern}" \
        --n-variants-suffix "${n_variants_suffix}" \
        --maf-cutoff ${maf_cutoff} --cis-window ${cis_window} \
        --fdr-threshold ${fdr_threshold} --pvalue-cutoff ${pvalue_cutoff} \
        --pvalue-col "${pvalue_col}" --af-col "${af_col}" \
        --study "${study}" --context "${context}" --genome "${genome}" \
        --output ${_output:a} --output-dir ${output_dir:a}