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(defaultcis_qtl[.]regional[.]tsv[.]gz$). Exampleoutput/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(defaultcis_qtl[.]pairs[.]tsv[.]gz$). Exampleoutput/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.752156227936135The recounted variant tallies, matched by
--n-variants-suffix(defaultcis_n_variants_stats[.]tsv[.]gz$). Exampleoutput/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: thecode/scriptdirectory holding the pecotmr_integration wrappers. Required.--maf-cutoffand--cis-window:0.01and1000000by 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 to0to skip the recount and use the counts as given.--pvalue-cutoff:0.05by 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.05by default, the threshold for the global significance call.--study,--context,--genome: labels stamped onto the assembledQtlSumStatsobject,study/context/hg38by default.--af-col,--pvalue-col,--molecular-id-col: the column names to read,af/pvalue/molecular_trait_object_idby 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 consolidatedQtlSumStatsobject holding the surviving xQTLs together with the--study,--contextand--genomelabels. This is the only file named in the workflow’soutput: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 hg38Command Interface¶
sos run pipeline/qtl_association_postprocessing.ipynb -husage: 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}
- Huang, Q. Q., Ritchie, S. C., Brozynska, M., & Inouye, M. (2018). Power, false discovery rate and Winner’s Curse in eQTL studies. Nucleic Acids Research, 46(22), e133–e133. 10.1093/nar/gky780