QTL Association Testing (TensorQTL) - FunGen-xQTL Consortium

Tests each molecular phenotype for association with nearby (cis) or genome-wide (trans) genetic variants using TensorQTL, reporting both nominal pair statistics and region-level significance.

Overview¶

The model is the standard QTL regression: each molecular phenotype is regressed on each candidate variant, with the covariate matrix (known covariates plus hidden factors) entering every fit. TensorQTL evaluates those regressions as batched tensor operations, on GPU when one is available and on CPU otherwise, which is what makes a genome-wide scan tractable. Two workflows expose that model.

cis tests each phenotype against the variants near it. The search region is a symmetric radius around the phenotype set by --window, 1 Mb by default, or an explicit per-region interval supplied through --customized-cis-windows, in which case --window is forced to zero. It runs as two chained steps: the first computes nominal statistics for every variant-phenotype pair on a chromosome and, unless --no-permutation is given, a permutation pass that fits a beta distribution to each region null; the second gathers the per-chromosome results into a single region-level table carrying permutation and beta-approximation p-values with q-values and Benjamini-Hochberg FDR.

trans tests phenotypes against variants far from them. A full scan is every phenotype against every variant, so memory rather than arithmetic is the binding constraint. The step therefore handles one phenotype chromosome at a time and loads a single genotype chromosome, named by --trans-geno-chromosome, so a genome-wide analysis is assembled from chromosome pairs rather than held in memory at once. Results are thinned as they are produced: --pval discards pairs above a nominal cutoff and records that cutoff in the output file name, while --pvalue-cutoff and --qvalue-cutoff set the significance filter applied to what remains. Permutation defaults to off here, because a null permuted within one chromosome is not the genome-wide null.

How rare a variant may be is controlled in both workflows by --MAC, a minor-allele-count floor, or by --maf-threshold, which overrides it with a frequency floor.

When to run it. After genotype, phenotype, and covariate preprocessing, which produce the by-chromosome file lists and the covariate matrix consumed here. The nominal summary statistics are the input to fine-mapping (mnm_regression) and to colocalization; the region-level table is what downstream enrichment and reporting read.

Input¶

Required:

--phenotype-file -- one bgzipped and tabixed molecular phenotype bed.gz, or a text file listing one such file per chromosome, e.g. output/phenotype/phenotype_by_chrom_for_cis/bulk_rnaseq.phenotype_by_chrom_files.txt.
--genotype-file -- a PLINK bed/bim/fam prefix, or a two-column list of chromosome ID and PLINK prefix, e.g. output/genotype_by_chrom/protocol_example.genotype.merged.plink_qc.genotype_by_chrom_files.txt. Only chromosomes present in both lists are analyzed, which is chr22 alone in the toy data.
--covariate-file -- the covariate matrix, known covariates plus hidden factors, e.g. output/covariate/protocol_example.rnaseq.bed.protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.Marchenko_PC.gz.
--name -- prefix for every output file, e.g. protocol_example.

Scope of the test:

--window -- cis radius in bp around each phenotype, default 1000000.
--customized-cis-windows -- four-column bed (chr, start, end, region ID) giving an explicit window per region; supplying it forces --window to zero.
--region-list and --region-list-phenotype-column -- restrict the analysis to the regions named in a file, reading names from the given column (default 4).
--trans-geno-chromosome -- for trans, the single genotype chromosome to load; without it the step uses the chromosome that came in with the phenotype.
--chromosome -- restrict to particular phenotype chromosomes.
--phenotype-group -- group file mapping molecular traits into trait objects, for cases such as sQTL where one region holds several events.
--keep-sample -- list of samples to retain.
--covariate-pattern -- keep only covariates whose names match these prefixes or exact names.
--interaction -- a file holding an interaction term when it is not already a covariate column; the name of its first column becomes the interaction name and is written into the output file names.

Thresholds:

--MAC -- minor-allele-count floor, default 0. --MAC 5 is used below because the toy cohort is small.
--maf-threshold -- minor-allele-frequency floor; overrides --MAC when set.
--permutation / --no-permutation -- region-level permutation testing, on by default for cis and off for trans.
--pval, --pval-threshold, --pvalue-cutoff, --qvalue-cutoff -- nominal and significance filters for trans.
--skip-nominal-if-exist -- reuse existing nominal parquet files instead of recomputing them.

Runtime:

--cwd -- output directory, default output.
--numThreads, --batch-size, --job-size, --walltime, --mem -- threads, variants per tensor batch, and cluster resources.
--container, --entrypoint -- software environment.
--modular-script-dir -- location of the Python driver, default code/script.

Example phenotype file¶

# Load required libraries
library(data.table)
library(genio)

# ------------------------------------------------------------------
# Example phenotype file (MWE: protocol_example, chr22)
# ------------------------------------------------------------------
# A text file listing one phenotype .bed.gz per chromosome.
pheno_path <- fread("output/phenotype/phenotype_by_chrom_for_cis/bulk_rnaseq.phenotype_by_chrom_files.txt")
head(pheno_path)

# One chromosome of phenotype data in BED-like format (.bed.gz).
pheno <- fread("output/phenotype/phenotype_by_chrom_for_cis/bulk_rnaseq.chr22.bed.gz")
pheno[1:5, 1:8]

Example genotype file¶

# ------------------------------------------------------------------
# Example genotype file (MWE: protocol_example)
# ------------------------------------------------------------------
# A text file with two columns: chromosome ID and PLINK prefix path.
geno_path <- fread(file.path("output/genotype_by_chrom/protocol_example.genotype.merged.plink_qc.genotype_by_chrom_files.txt"))
head(geno_path)

# Read the merged PLINK set (.bed/.bim/.fam) for the toy sample.
plink_data <- read_plink("output/plink/protocol_example.genotype.merged.plink_qc")
genotypes <- plink_data$X
genotypes[1:5, 1:5]

Reading: output/plink/protocol_example.genotype.merged.plink_qc.bim

Reading: output/plink/protocol_example.genotype.merged.plink_qc.fam

Reading: output/plink/protocol_example.genotype.merged.plink_qc.bed

Example covariates file¶

# ------------------------------------------------------------------
# Example covariates file (MWE: protocol_example)
# ------------------------------------------------------------------
cov <- fread("output/covariate/protocol_example.rnaseq.bed.protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.Marchenko_PC.gz")
cov[, 1:5]
dim(cov)

Example customized cis window file (Optional)¶

# ------------------------------------------------------------------
# Example customized cis-window file (Optional)
# ------------------------------------------------------------------
# Gene-level cis-window (start/end) defining the SNP search region per gene.
# Omit to use --window for a symmetric +/-1Mb window around each gene TSS.
window <- fread("tests/fixtures/generalized_TADB/expected/TADB_enhanced_cis.bed")
head(window)

Example interaction file¶

# ------------------------------------------------------------------
# Example interaction file (Optional, for iQTL analysis)
# ------------------------------------------------------------------
# Provide an interaction term as a separate file via --interaction
# when it is not already part of the covariate file.
int <- fread("input/ROSMAP_interaction_example.tsv")
head(int)

Error in fread("input/ROSMAP_interaction_example.tsv"): File 'input/ROSMAP_interaction_example.tsv' does not exist or is non-readable. getwd()=='/mnt/lustre/lab/gwang/home/rl3328/protocol_testing/xqtl-protocol_archived'
Traceback:

1. stopf("File '%s' does not exist or is non-readable. getwd()=='%s'", 
 .     file, getwd())
2. raise_condition(stop, gettextf(fmt, ..., domain = domain), c(class, 
 .     "simpleError", "error", "condition"))
3. signal(obj)

Output¶

<name>.<pheno_chr>_<geno_chr>.cis_qtl.pairs.tsv.gz (with .tbi) -- nominal statistics for every variant-phenotype pair tested in cis:

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

<name>.<pheno_chr>.cis_qtl_pairs.<chr>.parquet -- the same nominal results in the columnar form TensorQTL writes natively, kept so that a rerun with --skip-nominal-if-exist need not recompute them.
<name>.<pheno_chr>_<geno_chr>.cis_qtl_regional_significance.tsv.gz (with .tbi) -- one row per region carrying the top variant, the fitted beta parameters, permutation and beta-approximation p-values, q-values, FDR, and the per-region genomic inflation factor:
```
chrom	pos	n_variants	beta_shape1	beta_shape2	true_df	p_true_df	variant_id	tss_distance	tes_dista
22	11797697	6	0.9833031	5.6122484	30.858135	0.0676042590649028	chr22:11797697_C_T	858310	836360	
```
<name>.<pheno_chr>_<geno_chr>.cis_qtl_regional_significance.summary.txt -- how many traits were tested and how many pass each threshold, written by the collecting step:
```
TensorQTL cis-QTL regional significance summary
Input files:              1
Chromosomes:              1
Molecular traits tested:  200
Molecular trait objects:  200
```
<name>.<...>.cis_qtl.pairs.tsv.gz.cis_regional.<correction>.fdr.gz -- the region-level table repeated under each multiple-testing correction: permutation_adjusted, and Bonferroni with Benjamini-Hochberg over the original and the filtered variant sets.
<name>.<...>.maf_<maf>_window_<window>_cis_n_variants_stats.tsv.gz -- variants tested per region at the MAF and window actually applied:
```
chrom	molecular_trait_object_id	n_variants	n_variants_filtered
chr22	ENSG00000063515	44	44
chr22	ENSG00000075234	111	111
```
<name>.<pheno_chr>_geno_chr<chr>.trans_qtl.pairs.tsv.gz (with .tbi) -- trans statistics for the loaded genotype chromosome, with the phenotype and genotype chromosome recorded on each row; the file name gains _p_<cutoff> when --pval is set:
```
chrom	pos	variant_id	molecular_trait_id	pvalue	bhat	sebhat	r2	af	n	pheno_chrom	geno_chrom	qvalue
22	10685239	chr22:10685239_C_T	ENSG00000063515	0.202382946055137	-3.4470642	2.6448517	0.05358654
```
<name>.<pheno_chr>_geno_chr<chr>.trans_qtl.genomic_inflation.tsv.gz -- the genomic inflation factor of the trans scan for each molecular trait:
```
molecular_trait_id	genomic_inflation_lambda
ENSG00000063515	0.9919021540170271
ENSG00000075234	1.1570695153778894
```

Every step writes .stdout and .stderr beside its output. Example results for the toy data are under output/tensorqtl_cis/ and output/tensorqtl_trans/.

Result tables¶

For each chromosome, several summary statistics files are generated, including both nominal test statistics for each test and region (gene) level association evidence.

Nominal Association Results¶

The columns of the nominal association result are as follows:

chrom: Variant chromosome.
pos: Variant chromosomal position (basepairs).
molecular_trait_id: Molecular trait identifier (gene).
variant_id: ID of the variant (rsid or chr:position:ref:alt).
tss_distance: Distance of the SNP to the gene transcription start site (TSS).
tes_distance: Distance of the SNP to the gene transcription end site (TES).
cis_window_start_distance: Distance of the SNP to the start of the cis window (if using a customized cis window).
cis_window_end_distance: Distance of the SNP to the end of the cis window (if using a customized cis window).
af: The allele frequency of this SNP.
ma_samples: Number of samples carrying the minor allele.
ma_count: Total number of minor alleles across individuals.
pvalue: Nominal P-value from linear regression.
bhat: Slope of the linear regression.
sebhat: Standard error of bhat.
n: Number of phenotypes after basic QC.

Multiple Testing Corrected Results:¶

qvalue: Calculated q-value for each SNP (grouped by gene).

Interaction Association Results¶

Interaction results are produced when --interaction is set. The column names carry the name of the interaction variable, shown here for msex:

Model:

\text{phenotype} = \beta_0 + \beta_1 \cdot \text{snp} + \beta_2 \cdot \text{msex} + \beta_3 \cdot (\text{snp} \times \text{msex}) + \epsilon

(1)

(Taking msex as the interaction factor)

chrom: Chromosome number.
pos: Variant chromosomal position (basepairs).
a2: Variant reference allele (A, C, T, or G).
a1: Variant alternate allele.
molecular_trait_id: Molecular trait identifier, varies from phenotypes to phenotypes.
variant_id: ID of the top variant (rsid or chr:position:ref:alt).
af: Alternative allele frequency in the cohort analyzed.
ma_samples: Number of samples carrying the minor allele.
ma_count: Total number of minor alleles across individuals.
pvalue: P-value of the main effect from the interaction model.
bhat: Slope of the main effect from the interaction model.
se: Standard error of beta.
pvalue_msex: P-value of the msex term from the interaction model.
bhat_msex: Slope of the msex term from the interaction model.
se_msex: Standard error of bhat_msex.
pvalue_msex_interaction: P-value of the interaction term from the interaction model.
bhat_msex_interaction: Slope of the interaction term from the interaction model.
se_msex_interaction: Standard error of beta_msex_interaction.
molecular_trait_object_id: An intermediate ID (can be ignored).
n: Number of samples.

Multiple Testing Corrected Results:¶

qvalue_main: The q-value of the main effect.
qvalue_interaction: The q-value of the interaction effect.

Region (Gene) Level Association Evidence¶

The column specifications for region-level association evidence are as follows:

chrom: Chromosome number.
pos: Variant chromosomal position (basepairs).
n_variant: Total number of variants tested in cis.
beta_shape1: First parameter value of the fitted beta distribution.
beta_shape2: Second parameter value of the fitted beta distribution.
true_df: Effective degrees of freedom of the beta distribution approximation.
p_true_df: Empirical P-value for the beta distribution approximation.
variant_id: ID of the top variant (rsid or chr:position:ref:alt).
tss_distance: Distance of the SNP to the gene transcription start site (TSS).
tes_distance: Distance of the SNP to the gene transcription end site (TES).
ma_samples: Number of samples carrying the minor allele.
ma_count: Total number of minor alleles across individuals.
af: Alternative allele frequency.
p_nominal: Nominal P-value from linear regression.
bhat: Slope of the linear regression.
sehat: Standard error of the bhat.
p_perm: First permutation P-value directly obtained from the permutations with the direct method.
p_beta: Second permutation P-value obtained via beta approximation (this is the one to use for downstream analysis).
molecular_trait_object_id: Molecular trait identifier (gene).
n_traits: Group size in the permutation test.
genomic_inflation: Genomic inflation factor (lambda), quantifying the extent of bulk inflation and the excess false positive rate.

Multiple Testing Corrected Results:¶

q_beta: Q-value for p_beta using Storey’s method (qvalue), more conservative than FDR.
q_perm: Q-value for p_perm using Storey’s method (qvalue), more conservative than FDR.
fdr_beta: Adjusted P-value for p_beta using the Benjamini-Hochberg method (FDR).
fdr_perm: Adjusted P-value for p_perm using the Benjamini-Hochberg method (FDR).
p_nominal_threshold: Nominal p-value threshold for variants in the corresponding molecular trait, derived from empirical beta distribution as a result of permutation testing.

Minimal Working Example¶

The commands below run on the included toy data and write results under --cwd. The genotype, phenotype, and covariate inputs are produced by the preprocessing notebooks.

cis-QTL association¶

Test each molecular phenotype against variants within the cis window. Matching chromosomes between the genotype and phenotype lists are analyzed (chr22 in the toy data). --MAC 5 sets the minor-allele-count cutoff for the small toy sample.

Timing: TBD (on toy dataset)

sos run pipeline/TensorQTL.ipynb cis \
    --genotype-file output/genotype_by_chrom/protocol_example.genotype.merged.plink_qc.genotype_by_chrom_files.txt \
    --phenotype-file output/phenotype/phenotype_by_chrom_for_cis/bulk_rnaseq.phenotype_by_chrom_files.txt \
    --covariate-file output/covariate/protocol_example.rnaseq.bed.protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.Marchenko_PC.gz \
    --cwd output/tensorqtl_cis --name protocol_example --MAC 5 --numThreads 2

trans-QTL association¶

Test phenotypes against variants on a chosen genotype chromosome (--trans-geno-chromosome 22), restricted to the genes in --region-list. Trans analysis is memory-heavy genome-wide; here it is scoped to the toy chromosome.

Timing: TBD (on toy dataset)

sos run pipeline/TensorQTL.ipynb trans \
    --genotype-file output/genotype_by_chrom/protocol_example.genotype.merged.plink_qc.genotype_by_chrom_files.txt \
    --phenotype-file output/phenotype/phenotype_by_chrom_for_cis/bulk_rnaseq.phenotype_by_chrom_files.txt \
    --covariate-file output/covariate/protocol_example.rnaseq.bed.protocol_example.covariates.protocol_example.genotype.merged.plink_qc.plink_qc.prune.pca.Marchenko_PC.gz \
    --cwd output/tensorqtl_trans --name protocol_example --MAC 5 --numThreads 2 \
    --trans-geno-chromosome 22 --region-list data/combined_AD_genes.csv --region-list-phenotype-column 4

Command Interface¶

sos run pipeline/TensorQTL.ipynb -h

usage: sos run pipeline/TensorQTL.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:
  cis
  trans

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --cwd output (as path)
                        Path to the work directory of the analysis.
  --phenotype-file VAL (as path, required)
                        Phenotype file, or a list of phenotype per region.
  --genotype-file VAL (as path, required)
                        A genotype file in PLINK binary format (bed/bam/fam)
                        format, or a list of genotype per chrom
  --covariate-file VAL (as path, required)
                        Covariate file
  --covariate-pattern  (as list)
                        Optional pattern to filter covariates (list of covariate
                        prefixes or exact names)
  --name VAL (as str, required)
                        Prefix for the analysis output
  --region-list . (as path)
                        An optional subset of regions of molecular features to
                        analyze. The last column is the gene names
  --region-list-phenotype-column 4 (as int)
  --keep-sample . (as path)
                        Set list of sample to be keep
  --interaction ''
                        FIXME: please document
  --customized-cis-windows . (as path)
                        An optional list documenting the custom cis window for
                        each region to analyze, with four column, chr, start,
                        end, region ID (eg gene ID). If this list is not
                        provided, the default `window` parameter (see below)
                        will be used.
  --phenotype-group . (as path)
                        The phenotype group file to group molecule_trait into
                        molecule_trait_object This applies to multiple molecular
                        events in the same region, such as sQTL analysis.
  --chromosome  (as list)
                        The name of phenotype corresponding to gene_id or
                        gene_name in the region
  --MAC 0 (as int)
                        Minor allele count cutoff
  --trans-geno-chromosome ''
                        Specify genotype chromosome for trans analysis (e.g.
                        --trans-geno-chromosome 5) When set, trans step loads
                        this genotype chrom instead of the one from input_files
  --window 1000000 (as int)
                        Specify the cis window for the up and downstream radius
                        to analyze around the region of interest in units of bp
                        This parameter will be set to zero if
                        `customized_cis_windows` is provided.
  --numThreads 8 (as int)
                        Number of threads
  --job-size 1 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 12h
  --mem 16G
  --container ''
                        Container option for software to run the analysis:
                        docker or singularity
  --entrypoint ''
  --maf-threshold  MAC/(2.0*N)

                        Minor allele frequency cutoff. It will overwrite minor
                        allele cutoff. You may consider setting it to higher for
                        interaction analysis if you have statistical power
                        concerns
  --pvalue-cutoff '5e-8'
                        Filtering significant trans associations (for trans_2
                        workflow)
  --qvalue-cutoff ''

Sections
  cis_1:
    Workflow Options:
      --[no-]skip-nominal-if-exist (default to False)
                        parse input file lists skip nominal association results
                        if the files exists already This is false by default
                        which means to recompute everything This is only
                        relevant when the `parquet` files for nominal results
                        exist but not the other files and you want to avoid
                        computing the nominal results again
      --[no-]permutation (default to True)
  cis_2:
  trans:
    Workflow Options:
      --batch-size 10000 (as int)
      --pval-threshold 1.0 (as float)
      --[no-]permutation (default to False)
                        Permutation testing is incorrect when the analysis is
                        done by chrom
      --pval 0.0 (as float)

Workflow implementation¶

The SoS workflow definitions below are unchanged from the original protocol.

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
# Path to the work directory of the analysis.
parameter: cwd = path('output')
# Phenotype file, or a list of phenotype per region.
parameter: phenotype_file = path
# A genotype file in PLINK binary format (bed/bam/fam) format, or a list of genotype per chrom
parameter: genotype_file = path
# Covariate file
parameter: covariate_file = path
# Optional pattern to filter covariates (list of covariate prefixes or exact names)
parameter: covariate_pattern = []
# Prefix for the analysis output
parameter: name = str
# An optional subset of regions of molecular features to analyze. The last column is the gene names
parameter: region_list = path()
parameter: region_list_phenotype_column = 4
# Set list of sample to be keep
parameter: keep_sample = path()
# FIXME: please document
parameter: interaction = ""

# An optional list documenting the custom cis window for each region to analyze, with four column, chr, start, end, region ID (eg gene ID).
# If this list is not provided, the default `window` parameter (see below) will be used.
parameter: customized_cis_windows = path()

# The phenotype group file to group molecule_trait into molecule_trait_object
# This applies to multiple molecular events in the same region, such as sQTL analysis.
parameter: phenotype_group = path() 

# The name of phenotype corresponding to gene_id or gene_name in the region
parameter: chromosome = []
# Minor allele count cutoff
parameter: MAC = 0
# Specify genotype chromosome for trans analysis (e.g. --trans-geno-chromosome 5)
# When set, trans step loads this genotype chrom instead of the one from input_files
parameter: trans_geno_chromosome = ""

# Specify the cis window for the up and downstream radius to analyze around the region of interest in units of bp
# This parameter will be set to zero if `customized_cis_windows` is provided.
parameter: window = 1000000

# Number of threads
parameter: numThreads = 8
# For cluster jobs, number commands to run per job
parameter: job_size = 1
parameter: walltime = '12h'
parameter: mem = '16G'
# Container option for software to run the analysis: docker or singularity
parameter: container = ''
parameter: entrypoint = ''
import re

# Use the header of the covariate file to decide the sample size
import pandas as pd
N = len(pd.read_csv(covariate_file, sep = "\t",nrows = 1).columns) - 1

# Minor allele frequency cutoff. It will overwrite minor allele cutoff.
# You may consider setting it to higher for interaction analysis if you have statistical power concerns
parameter: maf_threshold = MAC/(2.0*N)

# Filtering significant trans associations (for trans_2 workflow)
parameter: pvalue_cutoff = "5e-8"
parameter: qvalue_cutoff = ""


import os
import pandas as pd

def adapt_file_path(file_path, reference_file):
    """
    Adapt a single file path based on its existence and a reference file's path.

    Args:
    - file_path (str): The file path to adapt.
    - reference_file (str): File path to use as a reference for adaptation.

    Returns:
    - str: Adapted file path.

    Raises:
    - FileNotFoundError: If no valid file path is found.
    """
    reference_path = os.path.dirname(reference_file)

    # Check if the file exists
    if os.path.isfile(file_path):
        return file_path

    # Check file name without path
    file_name = os.path.basename(file_path)
    if os.path.isfile(file_name):
        return file_name

    # Check file name in reference file's directory
    file_in_ref_dir = os.path.join(reference_path, file_name)
    if os.path.isfile(file_in_ref_dir):
        return file_in_ref_dir

    # Check original file path prefixed with reference file's directory
    file_prefixed = os.path.join(reference_path, file_path)
    if os.path.isfile(file_prefixed):
        return file_prefixed

    # If all checks fail, raise an error
    raise FileNotFoundError(f"No valid path found for file: {file_path}")

def adapt_file_path_all(df, column_name, reference_file):
    return df[column_name].apply(lambda x: adapt_file_path(x, reference_file))


if (str(genotype_file).endswith("bed") or str(genotype_file).endswith("pgen")) and str(phenotype_file).endswith("bed.gz"):
    input_files = [[phenotype_file, genotype_file]]
    if len(chromosome) > 0:
        input_chroms = [int(x) for x in chromosome]
    else:
        input_chroms = [0]
else:
    import pandas as pd
    import os
    molecular_pheno_files = pd.read_csv(phenotype_file, sep = "\t")
    
    if "#dir" in molecular_pheno_files.columns and "#chr" not in molecular_pheno_files.columns:
        molecular_pheno_files = molecular_pheno_files.rename(columns={"#dir": "path"})
        if "#id" in molecular_pheno_files.columns:
            molecular_pheno_files = molecular_pheno_files.rename(columns={"#id": "#chr"})
    
    if "#chr" in molecular_pheno_files.columns:
        molecular_pheno_files = molecular_pheno_files.groupby(['#chr','path']).size().reset_index(name='count').drop("count",axis = 1).rename(columns = {"#chr":"#id"})
    genotype_files = pd.read_csv(genotype_file,sep = "\t")
    genotype_files["#id"] = [x.replace("chr","") for x in genotype_files["#id"].astype(str)] # e.g. remove chr1 to 1
    genotype_files["#path"] = genotype_files["#path"].apply(lambda x: adapt_file_path(x, genotype_file))
    molecular_pheno_files["#id"] = [x.replace("chr","") for x in molecular_pheno_files["#id"].astype(str)]
    input_files = molecular_pheno_files.merge(genotype_files, on = "#id")
    
    # Only keep chromosome specified in --chromosome
    if len(chromosome) > 0:
        input_files = input_files[input_files['#id'].isin(chromosome)]
    input_files = input_files.values.tolist()
    input_chroms = [x[0] for x in input_files]
    input_files = [x[1:] for x in input_files]
    if len(name) == 0:
        name = f'{path(input_files[0][0]):bnn}' if len(input_files) == 1 else f'{path(input_files[0][0]):bnnn}'

[cis_1]
# parse input file lists
# skip nominal association results if the files exists already
# This is false by default which means to recompute everything
# This is only relevant when the `parquet` files for nominal results exist but not the other files
# and you want to avoid computing the nominal results again
parameter: skip_nominal_if_exist = False
parameter: permutation = True

# Extract interaction name
var_interaction = interaction
if os.path.isfile(interaction):
    interaction_s = pd.read_csv(interaction, sep='\t', index_col=0)
    var_interaction = interaction_s.columns[0] # interaction name
test_regional_association = permutation and len(var_interaction) == 0

input: input_files, group_by = len(input_files[0]), group_with = "input_chroms"
output_files = dict([("parquet", f'{cwd:a}/{_input[0]:bnn}{"_%s" % var_interaction if interaction else ""}.cis_qtl_pairs.{"" if input_chroms[_index] == 0 else input_chroms[_index]}.parquet'), # This convention is necessary to match the pattern of map_norminal output
                     ("nominal", f'{cwd:a}/{_input[0]:bnn}{"" if input_chroms[_index] == 0 else "_chr%s" % input_chroms[_index]}{"_%s" % var_interaction if interaction else ""}.cis_qtl.pairs.tsv.gz')])
if test_regional_association:
    output_files["regional"] = f'{cwd:a}/{_input[0]:bnn}{"" if input_chroms[_index] == 0 else "_chr%s" % input_chroms[_index]}{"_%s" % var_interaction if interaction else ""}.cis_qtl.regional.tsv.gz'
output: output_files
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output["nominal"]:bnnn}'
bash: expand= "${ }", stderr = f'{_output["nominal"]:nnn}.stderr', stdout = f'{_output["nominal"]:nnn}.stdout', container = container, entrypoint = entrypoint
    python3 ${modular_script_dir}/association_scan/TensorQTL/TensorQTL.py \
        --step cis \
        --cwd "${cwd}" \
        --genotype-file "${_input[1]}" \
        --phenotype-file "${_input[0]}" \
        --covariate-file "${covariate_file}" \
        --interaction "${interaction}" \
        --chromosome ${input_chroms[_index]} \
        --window ${window} \
        --MAC ${MAC} \
        --maf-threshold ${maf_threshold} \
        --pvalue-cutoff "${pvalue_cutoff}" \
        --qvalue-cutoff "${qvalue_cutoff}" \
        --permutation ${permutation} \
        ${"--skip-nominal-if-exist" if skip_nominal_if_exist else ""} \
        ${"--covariate-pattern " + " ".join([str(x) for x in covariate_pattern]) if len(covariate_pattern) else ""} \
        ${"--region-list " + str(region_list) + " --region-list-phenotype-column " + str(region_list_phenotype_column) if region_list.is_file() else ""} \
        ${"--keep-sample " + str(keep_sample) if keep_sample.is_file() else ""} \
        ${"--customized-cis-windows " + str(customized_cis_windows) if customized_cis_windows.is_file() else ""} \
        ${"--phenotype-group " + str(phenotype_group) if phenotype_group.is_file() else ""} \
        --numThreads ${numThreads}

    # Add per-gene Storey q-values to the raw nominal table (in place; re-bgzip+tabix)
    Rscript ${modular_script_dir}/association_scan/TensorQTL/TensorQTL.R \
        nominal_qvalues "${_output['nominal']}" "${interaction}"

[cis_2]
done_if("regional" not in _input.labels)
input: group_by = "all"
output_file_prefix = name if len(_input["nominal"]) > 1 else f'{_input["nominal"][0]:bnnnn}'
output: f'{cwd}/{output_file_prefix}.cis_qtl_regional_significance.tsv.gz',
        f'{cwd}/{output_file_prefix}.cis_qtl_regional_significance.summary.txt'
input_files = [str(x) for x in _input["regional"]]
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
bash: expand= "${ }", stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/association_scan/TensorQTL/TensorQTL.R \
        regional_postprocess "${_output[0]}" "${_output[1]}" ${" ".join(input_files)}

[trans]

parameter: batch_size = 10000
parameter: pval_threshold = 1.0
# Permutation testing is incorrect when the analysis is done by chrom
parameter: permutation = False
parameter: pval = 0.0

input: input_files, group_by = len(input_files[0]), group_with = "input_chroms"
output: nominal = f'{cwd:a}/{_input[0]:bnn}{"_chr%s" % input_chroms[_index] if input_chroms[_index] != 0 else ""}{"_geno_chr%s" % trans_geno_chromosome if trans_geno_chromosome else ""}.trans_qtl{"_p_%.0e" % pval if pval > 0.0 else ""}.pairs.tsv.gz'

trans_genotype_file = genotype_file if trans_geno_chromosome else _input[1]
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
bash: expand= "${ }", stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout', container = container, entrypoint = entrypoint
    python3 ${modular_script_dir}/association_scan/TensorQTL/TensorQTL.py \
        --step trans \
        --cwd "${cwd}" \
        --genotype-file "${trans_genotype_file}" \
        --phenotype-file "${_input[0]}" \
        --covariate-file "${covariate_file}" \
        --output "${_output['nominal']}" \
        --trans-geno-chromosome "${trans_geno_chromosome}" \
        --window ${window} \
        --MAC ${MAC} \
        --maf-threshold ${maf_threshold} \
        --batch-size ${batch_size} \
        --pval-threshold ${pval_threshold} \
        --pval ${pval} \
        --pvalue-cutoff "${pvalue_cutoff}" \
        --qvalue-cutoff "${qvalue_cutoff}" \
        ${"--covariate-pattern " + " ".join([str(x) for x in covariate_pattern]) if len(covariate_pattern) else ""} \
        ${"--region-list " + str(region_list) + " --region-list-phenotype-column " + str(region_list_phenotype_column) if region_list.is_file() else ""} \
        ${"--keep-sample " + str(keep_sample) if keep_sample.is_file() else ""} \
        ${"--customized-cis-windows " + str(customized_cis_windows) if customized_cis_windows.is_file() else ""} \
        --numThreads ${numThreads}

    # Add global Storey q-values to the raw trans table (in place; re-bgzip+tabix)
    Rscript ${modular_script_dir}/association_scan/TensorQTL/TensorQTL.R \
        trans_qvalues "${_output['nominal']}"

References¶

Taylor-Weiner, A., Aguet, F., Haradhvala, N. J., Gosai, S., Anand, S., Kim, J., Ardlie, K., Van Allen, E. M., & Getz, G. (2019). Scaling computational genomics to millions of individuals with GPUs. Genome Biology, 20(1). 10.1186/s13059-019-1836-7