vignettes/detectgxt.Rmd
detectgxt.RmdThe DetectGxT software provides functions for detecting gene-by-treatment (GxT) interactions on molecular count phenotypes, such as RNA-seq and ATAC-seq data. It aims to increase the statistical power to detect interaction (response) molecular QTLs using a nonlinear regression model that captures the relationship between genotypes and transformed molecular count phenotypes more accurately than a commonly used linear regression model.
This software is primarily designed for performing interaction QTL mapping focusing on a set of candidate pairs of features (e.g., gene expression and ATAC-peak regions) and SNPs for which significant genotype effects have been identified in at least one of the control and treatment conditions. Note that this type of multi-step approach has been adopted in previous studies, such as Alasoo et al. (2018) and Matoba et al. (2024).
A typical data analysis procedure using this package is as follows:
get_geno()
(optional)preprocess_pheno()
(optional)map_qtl()
Minimally, the input for the main function map_qtl() is
genotypes, pre-processed phenotypes, sequencing sample annotations, and
candidate feature-SNP pairs. See Input
data for details. Optionally, covariate information and a kinship
matrix may be supplied. The output is a data frame containing a row per
feature-SNP pair, with estimated coefficients and nominal
values for non-zero GxT interaction.
We first show the minimal code to run DetectGxT. See Input data for how to format the input.
# perform interaction QTL mapping
res <- map_qtl(
candidate=candidate, fn="nonlinear",
anno=anno, geno=geno, pheno=pheno)The minimal input for the main function map_qtl() is
data frames containing candidate feature-SNP pairs, genotypes, processed
molecular count phenotypes, and sample annotations. The formats of the
geno and pheno data frames are the same as
those for MatrixEQTL.
candidate: A data frame containing feature IDs and SNP
IDs in the first and second columns, respectively.geno: A data frame formatted as input for MatrixEQTL.
The first column must contain SNP IDs. The rest of the columns must
contain genotypes, which can be allelic dosage or integer-coded
genotypes (0, 1, 2). Each row corresponds to each SNP. Each column
corresponds to each sample.pheno: A data frame formatted as input for MatrixEQTL.
The first column must contain feature IDs (e.g., genes and ATAC-peak
regions). The rest of the columns must contain pre-processed molecular
phenotypes. Each row corresponds to each feature. Each column
corresponds to each sample.anno: A data frame containing sequencing sample IDs,
subject IDs, and treatment conditions (0 or 1) in the first, second, and
third columns, respectively. The columns must be named as “sample”,
“subject”, and “condition”.For the minimal code in Quick
start, we assume that the confounding factors have already been
regressed out of the phenotypes. We provide options to account for
confounding factors through the following input arguments to the
map_qtl() function.
covar: A data frame containing fixed effect covariates.
The first column must contain sample IDs. The rest of the columns must
contain covariate values. Each row corresponds to each sample. Each
column corresponds to each covariate.kinship: A matrix representing kinship (polygenic
structure). The row and column names of the matrix must correspond to
the set of unique elements in the subject column in
anno.The format of the covar data frame is the same as that
for MatrixEQTL.
In this vignette, we consider interaction (response) molecular QTL mapping on a pre-selected set of feature-SNP pairs obtained by condition-stratified molecular QTL mapping. Condition-stratified molecular QTL mapping can be performed based on linear regression or linear mixed-effect models using standard software packages, such as MatrixEQTL, TensorQTL, and limix_qtl. This step includes appropriate multiple testing correction and SNP filtering by selecting lead SNPs, performing LD clumping, or conducting conditional analysis. See Alasoo et al. (2018) and Matoba et al. (2024) for details.
Although assuming linearity between genotypes and transformed molecular count phenotypes can lead to inaccurate effect size estimation, this is unlikely to be a significant issue in detecting associations in a single condition. Therefore, we recommend using the standard linear regression approach in step 1 for computational efficiency.
It is necessary to obtain the union set of feature-SNP pairs from
step 1 and create a data frame that can be passed to
map_qtl() using the candidate input argument.
We do not provide a function to perform this task, as it depends on
upstream tools used in step 1. See Input
data for the format of the data frame.
In this section, we generate the anno data frame, which
is used throughout the remaining examples. See Input
data for the format of the data frame. For simplicity, we assume a
fully paired design, although the package also supports unpaired and
partially paired data.
n.sample <- 160 # number of sequencing samples
n.sub <- 80 # number of subjects
anno <- data.frame(
sample=paste0("s", seq_len(n.sample)),
subject=paste0("d", rep(seq_len(n.sub), each=2)),
condition=rep(c(0, 1), times=n.sub))
head(anno)
#> sample subject condition
#> 1 s1 d1 0
#> 2 s2 d1 1
#> 3 s3 d2 0
#> 4 s4 d2 1
#> 5 s5 d3 0
#> 6 s6 d3 1In this section, we demonstrate how to process genotype data in the
PLINK format to
generate a geno data frame, using get_geno(),
which can be passed to map_qtl() using the
geno argument. Users can pre-process the data using
alternative methods, if preferred.
The get_geno() function takes four arguments:
candidate: See Input
data.chromosomes: An integer vector of chromosome numbers to
search for the candidate SNPs.plink.prefix: The path prefix shared by all PLINK
files, up to and including the base.plink.suffix: The file extension of the PLINK binary
files (e.g., "bed").PLINK files must follow the naming convention
<plink.prefix>.<chromosome>.<plink.suffix>,
with .bim and .fam files sharing the same base
name. For example, with plink.prefix="path/to/geno",
chromosomes=c(1, 2), and plink.suffix="bed",
the function reads path/to/geno.1.bed and
path/to/geno.2.bed (along with the corresponding
.bim and .fam files). SNPs from all requested
chromosomes are stacked into a single output data frame.
The following example uses the PLINK files bundled with the package.
plink.prefix <- file.path(
system.file("extdata", "plink", package="detectgxt"), "geno")
candidate.geno <- data.frame(
feat_id=c("f1", "f2"),
snp_id=c("g1", "g3"))
geno.ex <- get_geno(
candidate=candidate.geno,
chromosomes=1,
plink.prefix=plink.prefix,
plink.suffix="bed")
head(geno.ex)
#> snp_id d1 d2 d3 d4 d5 d6
#> 1 g1 0 0 1 1 2 2
#> 2 g3 0 0 1 1 2 2In this section, we demonstrate how to pre-process raw molecular
count data to generate a pheno data frame, using
preprocess_pheno(), which can be passed to
map_qtl() using the pheno argument. Users can
pre-process the data using alternative methods, if preferred.
In the following section, we generate example count data
(count) using the make_count_data() function.
Users with existing molecular count data may skip to Input
to preprocess_pheno().
To create molecular count data, we first specify the number of features.
n.feat <- 10000 # number of features (e.g., genes, ATAC-peaks)We then use make_count_data(), which simulates
negative-binomial count data using the nonlinear GxT parametrization
consistent with the models in map_qtl(). By default, all
effects are set to zero via sd=c(0, 0, 0). To introduce
genotype, treatment, or interaction effects into the simulated counts,
specify sd (to draw per-feature effects from
Normal(0, sd[k]), where k indexes genotype,
treatment, and interaction) or coef (to fix effects at a
given value); exactly one of sd[k] or coef[k]
must be non-NA for each effect.
count <- make_count_data(
anno=anno, n.feat=n.feat, n.sample=n.sample)preprocess_pheno()
The following code displays the count data generated in the previous section.
count[1:5, 1:5]
#> s1 s2 s3 s4 s5
#> f1 4 3 12 0 15
#> f2 60 126 173 76 63
#> f3 0 0 0 0 0
#> f4 57 106 66 32 76
#> f5 365 253 225 119 143preprocess_pheno()
The preprocess_pheno() function first scales and
transforms the data. Then it residualizes the data against the treatment
indicators, performs principal component analysis, and residualizes the
transformed data against a desired number of principal components.
The preprocess_pheno() function takes arguments
including:
count: A matrix containing molecular counts
(non-negative integers). The row and column names correspond to feature
IDs and sample IDs, respectively.anno: See Input
data.num.pc: An integer specifying the number of molecular
phenotype principal components to control for.The following code runs the function using the anno and
count objects generated above.
pheno <- preprocess_pheno(
count=count, anno=anno, num.pc=10)
pheno[1:5, 1:5]
#> feat_id s1 s2 s3 s4
#> 1 f1 0.52375774 0.1677593 0.5179037 -1.43289457
#> 2 f2 -0.41880813 0.1979051 0.6057498 -0.03943474
#> 3 f3 -0.33996767 -0.2226068 -0.2359433 -0.39787193
#> 4 f4 0.07208729 0.7853798 0.2156213 -0.37325833
#> 5 f5 0.72546085 0.4999999 0.1621861 -0.41803018Note that the pheno object above is only for
demonstrating the usage of preprocess_pheno() and will not
be used in the subsequent section, as it does not include genotype or
treatment effects.
The number of principal components can be adjusted as described in
Matoba et al. (2024). Alternatively, a
scree plot or an automatic elbow detection method can be used to select
an appropriate number. In the latter case, the list
argument can be set to TRUE to obtain the output from the
prcomp() function from the stats package.
pheno <- preprocess_pheno(
count=count, anno=anno, num.pc=10, list=TRUE)
names(pheno)
#> [1] "pheno" "pca"In this case, the output is a list object containing the
pheno data frame and the pca object returned
from prcomp().
map_qtl()
In this section, we demonstrate how to use map_qtl().
The function takes candidate, geno,
pheno, and anno as required input and returns
a data frame containing effect sizes and
values.
In the following section, we generate example input data
(candidate, geno, pheno) using
the make_sim_data() function. Users with existing
pre-processed phenotype data may skip to Input
to map_qtl().
make_sim_data()
In this section, we generate genotype and phenotype data using the
make_sim_data() function, which will be used for
demonstrating the usage of the map_qtl() function. Note
that make_sim_data() returns simulated phenotype data that
mimics preprocessed data.
We will draw genotype, treatment, and interaction effect sizes from Normal distributions with zero mean and effect-specific standard deviations. The following code specifies the standard deviations.
sd <- c(1.5, 2.0, 1.0) # genotype, treatment, interactionWe will generate data based on eight models corresponding to all combinations of the genotype, treatment, and interaction effects being either zero or non-zero. Specifically, we define models 1-8 as follows.
| Model | Genotype | Treatment | Interaction |
|---|---|---|---|
| 1 | 0 | 0 | 0 |
| 2 | 1 | 0 | 0 |
| 3 | 0 | 1 | 0 |
| 4 | 1 | 1 | 0 |
| 5 | 0 | 0 | 1 |
| 6 | 1 | 0 | 1 |
| 7 | 0 | 1 | 1 |
| 8 | 1 | 1 | 1 |
Here, 0 and 1 indicate that the effect is zero and non-zero,
respectively. The number of replicates per model can be specified using
the num input argument, which expects an eight-dimensional
integer vector. With the following vector, 900 simulations will be based
on the fourth model, and 100 simulations will be based on the eighth
model.
num <- c(0, 0, 0, 900, 0, 0, 0, 100)The following code generates data for 1000 feature-SNP pairs.
d <- make_sim_data(
num=num, anno=anno, fn="nonlinear", sd=sd, filter.geno=TRUE)Note that, if filter.geno is set to TRUE,
SNPs for which either of the homozygotes is missing in at least one of
the control and treated conditions are excluded. See Running
map_qtl() for details.
map_qtl()
The following code displays the input data generated in the previous
section. Note that only the first two columns are required in the
candidate data frame, while this example contains
additional columns.
head(d$candidate) # candidate
#> feat_id snp_id model maf b0 b1 b2 b3
#> 1 f1 g1 4 0.1694789 0 0.8411191 -0.9055679 0
#> 2 f2 g2 4 0.3756608 0 -1.3180614 2.3259291 0
#> 3 f3 g3 4 0.2091136 0 -2.2617105 -4.4458923 0
#> 4 f4 g4 4 0.4745412 0 -0.4854759 -0.5104962 0
#> 5 f5 g5 4 0.3371473 0 0.7503507 -2.3098373 0
#> 6 f6 g6 4 0.3786667 0 3.6018332 1.6039236 0
d$geno[1:5, 1:5] # geno
#> snp_id d1 d2 d3 d4
#> 1 g1 0 0 1 0
#> 2 g2 0 1 0 2
#> 3 g3 0 0 2 0
#> 4 g4 2 0 1 0
#> 5 g5 1 1 0 1
d$pheno[1:5, 1:5] # pheno
#> feat_id s1 s2 s3 s4
#> 1 f1 -0.83204330 -2.0721385 -1.0655906 -2.469350
#> 2 f2 -2.00016494 1.7811384 -0.8796291 1.535850
#> 3 f3 -1.17867244 -6.2286525 -0.9879191 -3.717404
#> 4 f4 -0.91003067 -2.3049397 1.8297305 -1.940412
#> 5 f5 0.08363835 -0.9590932 0.4192116 -2.129436map_qtl()
We now perform interaction QTL mapping using map_qtl().
See Required
for the required input arguments, candidate,
geno, pheno, and anno. As the
phenotype is assumed to have been residualized against covariates, we
omit the covar input argument (see Optional).
We omit kinship, as we do not include random effects in
this example. See Including
random effects for how to include random effects.
Additional arguments include:
fn: A character string specifying the function. This
must be one of “nonlinear” and “linear”, corresponding to nonlinear and
linear models, respectively.filter.geno: A Boolean variable as to whether to ensure
that homozygotes for both alleles have at least one observation in both
control and treated conditions.
res <- map_qtl(
candidate=d$candidate, fn="nonlinear",
geno=d$geno, pheno=d$pheno, anno=anno,
filter.geno=TRUE)
head(res)
#> feat_id snp_id beta0 beta_g beta_t beta_gxt sigma sigma_u
#> 1 f1 g1 -0.08506611 0.9250519 -1.0026361 0.1488315 0.9223795 NA
#> 2 f2 g2 0.15370802 -1.3154423 2.1059396 0.0680372 0.9793602 NA
#> 3 f3 g3 -0.20450228 -2.4422318 -4.1314013 -0.1130397 0.9585258 NA
#> 4 f4 g4 -0.16305591 -0.3926680 -0.4118585 -0.1417504 1.0866933 NA
#> 5 f5 g5 -0.36569753 0.9718277 -1.9568060 -0.2670066 1.0888461 NA
#> 6 f6 g6 -0.21291378 3.6086383 1.7850510 0.0136696 1.0329443 NA
#> p_value return_code
#> 1 0.3968035 0
#> 2 0.7553684 0
#> 3 0.6964196 0
#> 4 0.5399299 0
#> 5 0.1748733 0
#> 6 0.9346815 0The output is a data frame containing feature IDs, SNP IDs, effect
estimates, and
values for non-zero GxT interaction. Note that sigma_u,
which represents estimates of random effect standard deviations, is
NA, as we did not include random effects in this
example.
We emphasize that we set filter.geno to
TRUE to exclude SNPs for which either of the homozygotes is
missing in at least one of the control and treated conditions, which we
call the “AA-AB” SNPs. We generally recommend excluding these SNPs, as
it is often impossible to reliably estimate genotype effects with
nonlinear regression. See our
manuscript for further details. In practice, there is no need to
pass this argument explicitly, as this is the default setting. Users can
filter data differently before running map_qtl(), if
preferred. In this case, filter.geno must be set to
FALSE.
In this section, we demonstrate how to fit a nonlinear mixed effect model including polygenic or subject random effects.
The following code generates an example kinship matrix. We use an
identity matrix, which controls for the subject random effect, for
simplicity. The row and column names of the matrix must correspond to
the set of unique elements in the subject column in
anno.
We then run map_qtl(), specifying the input arguments
ranef and kinship. Note that the
kinship argument can be omitted in practice when the matrix
is identity.
res.ranef <- map_qtl(
candidate=d$candidate, fn="nonlinear",
geno=d$geno, pheno=d$pheno, anno=anno,
ranef=TRUE, kinship=kinship, filter.geno=TRUE)
head(res.ranef)
#> feat_id snp_id beta0 beta_g beta_t beta_gxt sigma
#> 1 f1 g1 -0.08506785 0.9250610 -1.0026591 0.14887202 0.9069409
#> 2 f2 g2 0.15370807 -1.3154424 2.1059396 0.06803719 0.9793291
#> 3 f3 g3 -0.20450225 -2.4422319 -4.1314014 -0.11303948 0.9584988
#> 4 f4 g4 -0.16409346 -0.3912414 -0.4103359 -0.14390711 1.0497373
#> 5 f5 g5 -0.36531241 0.9714660 -1.9564032 -0.26743082 1.0642698
#> 6 f6 g6 -0.21291380 3.6086383 1.7850510 0.01366960 1.0329149
#> sigma_u p_value return_code
#> 1 0.168088458 0.3893308 0
#> 2 0.008207871 0.7553678 0
#> 3 0.007394676 0.6964162 0
#> 4 0.280995849 0.5230238 0
#> 5 0.230034960 0.1662258 0
#> 6 0.008340928 0.9346796 0The output now contains values in the sigma_u column,
which represents estimates of random effect standard deviations.
sessionInfo()
#> R version 4.4.3 (2025-02-28)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.6.1
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] detectgxt_0.1.1
#>
#> loaded via a namespace (and not attached):
#> [1] Matrix_1.7-4 jsonlite_2.0.0 dplyr_1.2.0
#> [4] compiler_4.4.3 tidyselect_1.2.1 jquerylib_0.1.4
#> [7] splines_4.4.3 systemfonts_1.3.1 textshaping_1.0.4
#> [10] yaml_2.3.12 fastmap_1.2.0 lattice_0.22-9
#> [13] R6_2.6.1 generics_0.1.4 knitr_1.51
#> [16] BiocGenerics_0.52.0 htmlwidgets_1.6.4 tibble_3.3.1
#> [19] desc_1.4.3 snpStats_1.56.0 bslib_0.10.0
#> [22] pillar_1.11.1 rlang_1.1.7 cachem_1.1.0
#> [25] xfun_0.56 fs_1.6.6 sass_0.4.10
#> [28] otel_0.2.0 cli_3.6.5 withr_3.0.2
#> [31] pkgdown_2.2.0 magrittr_2.0.4 zlibbioc_1.52.0
#> [34] digest_0.6.39 grid_4.4.3 rstudioapi_0.18.0
#> [37] mvtnorm_1.3-6 lifecycle_1.0.5 vctrs_0.7.2
#> [40] evaluate_1.0.5 glue_1.8.0 ragg_1.5.0
#> [43] survival_3.8-6 rmarkdown_2.30 purrr_1.2.1
#> [46] tools_4.4.3 pkgconfig_2.0.3 htmltools_0.5.9