Introduction to CSeQTL
Paul Little
Last Refresh: 2025-03-04
Overview
Assuming all software dependencies and CSeQTL (Little et al., 2023) installation are installed, we can begin.
# Load libraries
req_packs = c("ggplot2","smarter","CSeQTL")
for(pack in req_packs){
library(package = pack,character.only = TRUE)
}
# List package's exported functions
ls("package:CSeQTL")
#> [1] "CSeQTL_GS" "CSeQTL_dataGen" "CSeQTL_full_analysis"
#> [4] "CSeQTL_linearTest" "CSeQTL_oneExtremeSim" "CSeQTL_run_MatrixEQTL"
#> [7] "CSeQTL_smart" "OLS_sim" "gen_true_RHO"
#> [10] "plot_RHO" "prep_gene_info"
# Fix seed
set.seed(2)Simulation: No eQTL
We will simulate a dataset with three cell types, reference allele differential expression, and no cell type-specific eQTLs.
# sample size
NN = 3e2
# fold-change between q-th and 1st cell type
true_KAPPA = c(1,3,1)
# eQTL effect size per cell type,
# fold change between B and A allele
true_ETA = c(1,1,1)
# cis/trans effect size
true_ALPHA = c(1,1,1)
# count number of cell types
QQ = length(true_KAPPA)
# TReC model overdispersion
true_PHI = 0.1
# ASReC model overdispersion
true_PSI = 0.05
# Simulate cell type proportions
tRHO = gen_true_RHO(wRHO = 1,NN = NN,QQ = QQ)
plot_RHO(RHO = tRHO)
# Simulate a data object for a gene and SNP
sim = CSeQTL_dataGen(NN = NN,MAF = 0.3,true_BETA0 = log(1000),
true_KAPPA = true_KAPPA,true_ETA = true_ETA,true_PHI = true_PHI,
true_PSI = true_PSI,prob_phased = 0.05,true_ALPHA = true_ALPHA,
RHO = tRHO,cnfSNP = TRUE,show = FALSE)
names(sim)
#> [1] "true_PARS" "true_SNP" "dat" "XX" "true_RHO" "QQ"
#> [7] "GI" "np" "I_np" "iBETA" "iPHI" "MU"
#> [13] "true_OF" "vname"
# TReC, ASReC, haplotype 2 counts
sim$dat[1:3,]
#> total LGX1 total_phased hap2 LBC phased log_mu mu pp
#> 1 456 2339.837 27 13 16.81415 1 6.209945 497.6740 0.5
#> 2 1080 6467.905 47 17 28.63941 1 6.973643 1068.1067 0.5
#> 3 702 3903.057 30 10 17.21821 1 6.566389 710.7982 0.5
#> phase0 log_lib_size geno_col
#> 1 0 5.538577 #FF0000BF
#> 2 0 5.850341 #0000FFBF
#> 3 0 5.751526 #00FF00BF
# Batch covariates including the intercept
sim$XX[1:3,]
#> X1 X2 X3 X4 X5
#> [1,] 1 -1.7890676 1 -0.5523594 1.35690857
#> [2,] 1 -0.6066932 0 1.5950620 1.06844120
#> [3,] 1 -0.9814495 0 -0.9546107 0.06637533
sim$dat$SNP = sim$true_SNP
sim$dat$SNP = factor(sim$dat$SNP,
levels = sort(unique(sim$dat$SNP)),
labels = c("AA","AB","BA","BB"))
# TReC vs SNP
ggplot(data = sim$dat,
mapping = aes(x = SNP,y = log10(total + 1))) +
geom_violin(aes(fill = SNP)) + geom_jitter(width = 0.25) +
geom_boxplot(width = 0.1) +
xlab("Phased Genotype") + ylab("log10(TReC + 1)") +
theme(legend.position = "right",
axis.title = element_text(size = 15),
axis.text = element_text(size = 12))
# TReC vs ASReC
ggplot(data = sim$dat,
mapping = aes(x = total,y = total_phased,color = SNP)) +
geom_point() + xlab("Total Read Count (TReC)") +
ylab("Total Allele-specific Read Count (ASReC)") +
theme(legend.position = "right",
axis.title = element_text(size = 15),
axis.text = element_text(size = 12))
CSeQTL
Bulk Test
Based on the above TReC boxplot, the phased SNP genotype appears to
be an eQTL. However, this has yet to account for batch effects. Let us
test for a bulk eQTL accounting for batch covariates
sim$XX. Make sure to center continuous covariates. Notice
that CSeQTL, treating bulk TReC as sourced from a single cell type,
corresponds to the TReCASE (Sun,
2012) method.
summary(sim$XX)
#> X1 X2 X3 X4 X5
#> Min. :1 Min. :-2.90650 Min. :0.0000 Min. :-1.75525 Min. :-2.48349
#> 1st Qu.:1 1st Qu.:-0.64850 1st Qu.:0.0000 1st Qu.:-0.83890 1st Qu.:-0.65859
#> Median :1 Median :-0.02574 Median :0.0000 Median : 0.05764 Median :-0.00457
#> Mean :1 Mean : 0.00000 Mean :0.4833 Mean : 0.00000 Mean : 0.00000
#> 3rd Qu.:1 3rd Qu.: 0.68804 3rd Qu.:1.0000 3rd Qu.: 0.85752 3rd Qu.: 0.69906
#> Max. :1 Max. : 2.82172 Max. :1.0000 Max. : 1.67543 Max. : 2.69418
mRHO = matrix(1,NN,1)
colnames(mRHO) = "Bulk"
cistrans_thres = 0.01
res = c()
for(trec_only in c(TRUE,FALSE)){
for(neg_binom in c(TRUE,FALSE)){
for(beta_binom in c(TRUE,FALSE)){
cat(sprintf("%s: trec_only = %s, neg_binom = %s, beta_binom = %s ...\n",
date(),trec_only,neg_binom,beta_binom))
PHASE = sim$dat$phased * ifelse(trec_only,0,1)
upPHI = ifelse(neg_binom,1,0)
upPSI = ifelse(beta_binom,1,0)
sout = CSeQTL_smart(TREC = sim$dat$total,hap2 = sim$dat$hap2,
ASREC = sim$dat$total_phased,PHASE = PHASE,
SNP = sim$true_SNP,RHO = mRHO,XX = sim$XX,upPHI = upPHI,
upKAPPA = 1,upETA = 1,upPSI = upPSI,upALPHA = 1,
iFullModel = FALSE,trim = FALSE,thres_TRIM = 10,
hypotest = TRUE,swap = FALSE,numAS = 5,numASn = 5,
numAS_het = 5,cistrans_thres = cistrans_thres)
# sout$HYPO
res = rbind(res,smart_df(Model = ifelse(trec_only,"TReC-only","TReCASE"),
TReC_Dist = ifelse(neg_binom,"Negative Binomial","Poisson"),
ASReC_Dist = ifelse(beta_binom,"Beta-Binomial","Binomial"),
sout$res))
rm(sout)
}}}
#> Tue Mar 4 17:10:55 2025: trec_only = TRUE, neg_binom = TRUE, beta_binom = TRUE ...
#> Tue Mar 4 17:10:55 2025: trec_only = TRUE, neg_binom = TRUE, beta_binom = FALSE ...
#> Tue Mar 4 17:10:55 2025: trec_only = TRUE, neg_binom = FALSE, beta_binom = TRUE ...
#> Tue Mar 4 17:10:55 2025: trec_only = TRUE, neg_binom = FALSE, beta_binom = FALSE ...
#> Tue Mar 4 17:10:55 2025: trec_only = FALSE, neg_binom = TRUE, beta_binom = TRUE ...
#> Tue Mar 4 17:10:56 2025: trec_only = FALSE, neg_binom = TRUE, beta_binom = FALSE ...
#> Tue Mar 4 17:10:56 2025: trec_only = FALSE, neg_binom = FALSE, beta_binom = TRUE ...
#> Tue Mar 4 17:10:56 2025: trec_only = FALSE, neg_binom = FALSE, beta_binom = FALSE ...
num_vars = which(unlist(lapply(res,function(xx) class(xx))) == "numeric")
res[,num_vars] = apply(res[,num_vars],2,function(xx) smart_SN(x = xx,digits = 2))
res$Interpret = apply(res[,c("cistrans","PVAL_eQTL")],1,function(xx){
ct = xx[1]
pval = as.numeric(xx[2])
ifelse(pval < 0.05,sprintf("%s eQTL",ct),"no eQTL")
})
res$Correct_Model = apply(res[,c("TReC_Dist","ASReC_Dist")],1,function(xx){
chk_trec = (( true_PHI > 0 & xx[1] == "Negative Binomial" )
| ( true_PHI == 0 & xx[1] == "Poisson" ))
chk_trec
chk_asrec = (( true_PSI > 0 & xx[2] == "Beta-Binomial" )
| ( true_PSI == 0 & xx[2] == "Binomial" ))
chk_trec = ifelse(chk_trec,"TReC-Yes","TReC-No")
chk_asrec = ifelse(chk_asrec,"ASReC-Yes","ASReC-No")
sprintf("%s;%s",chk_trec,chk_asrec)
})
# Simplified output
smart_rmcols(res,c("LRT_trec","LRT_trecase","LRT_cistrans","cistrans"))
#> Model TReC_Dist ASReC_Dist CT PVAL_trec PVAL_trecase PVAL_cistrans PVAL_eQTL Interpret Correct_Model
#> 1 TReC-only Negative Binomial Beta-Binomial Bulk 0.00e+00 0.00e+00 0.00e+00 0.00e+00 trans eQTL TReC-Yes;ASReC-Yes
#> 2 TReC-only Negative Binomial Binomial Bulk 0.00e+00 0.00e+00 0.00e+00 0.00e+00 trans eQTL TReC-Yes;ASReC-No
#> 3 TReC-only Poisson Beta-Binomial Bulk 0.00e+00 0.00e+00 0.00e+00 0.00e+00 trans eQTL TReC-No;ASReC-Yes
#> 4 TReC-only Poisson Binomial Bulk 0.00e+00 0.00e+00 0.00e+00 0.00e+00 trans eQTL TReC-No;ASReC-No
#> 5 TReCASE Negative Binomial Beta-Binomial Bulk 0.00e+00 0.00e+00 0.00e+00 0.00e+00 trans eQTL TReC-Yes;ASReC-Yes
#> 6 TReCASE Negative Binomial Binomial Bulk 0.00e+00 2.80e-04 0.00e+00 0.00e+00 trans eQTL TReC-Yes;ASReC-No
#> 7 TReCASE Poisson Beta-Binomial Bulk 0.00e+00 0.00e+00 0.00e+00 0.00e+00 trans eQTL TReC-No;ASReC-Yes
#> 8 TReCASE Poisson Binomial Bulk 0.00e+00 0.00e+00 0.00e+00 0.00e+00 trans eQTL TReC-No;ASReC-NoThus a bulk approach to eQTL testing without accounting for cell type variation leads to a false positive association. Model misspecification can also lead to inflated type 1 error.
Cell type-specific Testing
The code below adjusts for cell type proportions and tests for cell type-specific eQTLs.
cistrans_thres = 0.01
res = c()
for(trec_only in c(TRUE,FALSE)){
for(neg_binom in c(TRUE,FALSE)){
for(beta_binom in c(TRUE,FALSE)){
cat(sprintf("%s: trec_only = %s, neg_binom = %s, beta_binom = %s ...\n",
date(),trec_only,neg_binom,beta_binom))
PHASE = sim$dat$phased * ifelse(trec_only,0,1)
upPHI = ifelse(neg_binom,1,0)
upPSI = ifelse(beta_binom,1,0)
upKAPPA = rep(1,QQ)
upETA = upKAPPA
upALPHA = upETA
sout = CSeQTL_smart(TREC = sim$dat$total,hap2 = sim$dat$hap2,
ASREC = sim$dat$total_phased,PHASE = PHASE,
SNP = sim$true_SNP,RHO = sim$true_RHO,XX = sim$XX,upPHI = upPHI,
upKAPPA = upKAPPA,upETA = upETA,upPSI = upPSI,upALPHA = upALPHA,
iFullModel = FALSE,trim = FALSE,thres_TRIM = 10,
hypotest = TRUE,swap = FALSE,numAS = 5,numASn = 5,
numAS_het = 5,cistrans_thres = cistrans_thres)
# sout$HYPO
res = rbind(res,smart_df(Model = ifelse(trec_only,"TReC-only","TReCASE"),
TReC_Dist = ifelse(neg_binom,"Negative Binomial","Poisson"),
ASReC_Dist = ifelse(beta_binom,"Beta-Binomial","Binomial"),
sout$res))
rm(sout)
}}}
#> Tue Mar 4 17:10:57 2025: trec_only = TRUE, neg_binom = TRUE, beta_binom = TRUE ...
#> Tue Mar 4 17:10:57 2025: trec_only = TRUE, neg_binom = TRUE, beta_binom = FALSE ...
#> Tue Mar 4 17:10:57 2025: trec_only = TRUE, neg_binom = FALSE, beta_binom = TRUE ...
#> Tue Mar 4 17:10:57 2025: trec_only = TRUE, neg_binom = FALSE, beta_binom = FALSE ...
#> Tue Mar 4 17:10:57 2025: trec_only = FALSE, neg_binom = TRUE, beta_binom = TRUE ...
#> Tue Mar 4 17:10:58 2025: trec_only = FALSE, neg_binom = TRUE, beta_binom = FALSE ...
#> Tue Mar 4 17:10:59 2025: trec_only = FALSE, neg_binom = FALSE, beta_binom = TRUE ...
#> Tue Mar 4 17:11:00 2025: trec_only = FALSE, neg_binom = FALSE, beta_binom = FALSE ...
num_vars = which(unlist(lapply(res,function(xx) class(xx))) == "numeric")
res[,num_vars] = apply(res[,num_vars],2,function(xx) smart_SN(x = xx,digits = 2))
res$Interpret = apply(res[,c("cistrans","PVAL_eQTL")],1,function(xx){
ct = xx[1]
pval = as.numeric(xx[2])
ifelse(pval < 0.05,sprintf("%s eQTL",ct),"no eQTL")
})
res$Correct_Model = apply(res[,c("TReC_Dist","ASReC_Dist")],1,function(xx){
chk_trec = (( true_PHI > 0 & xx[1] == "Negative Binomial" )
| ( true_PHI == 0 & xx[1] == "Poisson" ))
chk_trec
chk_asrec = (( true_PSI > 0 & xx[2] == "Beta-Binomial" )
| ( true_PSI == 0 & xx[2] == "Binomial" ))
chk_trec = ifelse(chk_trec,"TReC-Yes","TReC-No")
chk_asrec = ifelse(chk_asrec,"ASReC-Yes","ASReC-No")
sprintf("%s;%s",chk_trec,chk_asrec)
})
# Simplified output
smart_rmcols(res,c("LRT_trec","LRT_trecase","LRT_cistrans","cistrans"))
#> Model TReC_Dist ASReC_Dist CT PVAL_trec PVAL_trecase PVAL_cistrans PVAL_eQTL Interpret Correct_Model
#> 1 TReC-only Negative Binomial Beta-Binomial CT1 9.12e-01 9.12e-01 0.00e+00 9.12e-01 no eQTL TReC-Yes;ASReC-Yes
#> 2 TReC-only Negative Binomial Beta-Binomial CT2 7.06e-01 7.06e-01 0.00e+00 7.06e-01 no eQTL TReC-Yes;ASReC-Yes
#> 3 TReC-only Negative Binomial Beta-Binomial CT3 2.96e-01 2.96e-01 0.00e+00 2.96e-01 no eQTL TReC-Yes;ASReC-Yes
#> 4 TReC-only Negative Binomial Binomial CT1 9.12e-01 9.12e-01 0.00e+00 9.12e-01 no eQTL TReC-Yes;ASReC-No
#> 5 TReC-only Negative Binomial Binomial CT2 7.06e-01 7.06e-01 0.00e+00 7.06e-01 no eQTL TReC-Yes;ASReC-No
#> 6 TReC-only Negative Binomial Binomial CT3 2.96e-01 2.96e-01 0.00e+00 2.96e-01 no eQTL TReC-Yes;ASReC-No
#> 7 TReC-only Poisson Beta-Binomial CT1 7.07e-01 7.07e-01 0.00e+00 7.07e-01 no eQTL TReC-No;ASReC-Yes
#> 8 TReC-only Poisson Beta-Binomial CT2 0.00e+00 0.00e+00 0.00e+00 0.00e+00 trans eQTL TReC-No;ASReC-Yes
#> 9 TReC-only Poisson Beta-Binomial CT3 0.00e+00 0.00e+00 0.00e+00 0.00e+00 trans eQTL TReC-No;ASReC-Yes
#> 10 TReC-only Poisson Binomial CT1 7.07e-01 7.07e-01 0.00e+00 7.07e-01 no eQTL TReC-No;ASReC-No
#> 11 TReC-only Poisson Binomial CT2 0.00e+00 0.00e+00 0.00e+00 0.00e+00 trans eQTL TReC-No;ASReC-No
#> 12 TReC-only Poisson Binomial CT3 0.00e+00 0.00e+00 0.00e+00 0.00e+00 trans eQTL TReC-No;ASReC-No
#> 13 TReCASE Negative Binomial Beta-Binomial CT1 9.12e-01 7.77e-01 9.11e-01 7.77e-01 no eQTL TReC-Yes;ASReC-Yes
#> 14 TReCASE Negative Binomial Beta-Binomial CT2 7.06e-01 7.68e-01 8.04e-01 7.68e-01 no eQTL TReC-Yes;ASReC-Yes
#> 15 TReCASE Negative Binomial Beta-Binomial CT3 2.96e-01 3.73e-01 5.39e-01 3.73e-01 no eQTL TReC-Yes;ASReC-Yes
#> 16 TReCASE Negative Binomial Binomial CT1 9.12e-01 6.04e-01 9.60e-01 6.04e-01 no eQTL TReC-Yes;ASReC-No
#> 17 TReCASE Negative Binomial Binomial CT2 7.06e-01 6.84e-01 6.08e-01 6.84e-01 no eQTL TReC-Yes;ASReC-No
#> 18 TReCASE Negative Binomial Binomial CT3 2.96e-01 1.91e-01 5.32e-01 1.91e-01 no eQTL TReC-Yes;ASReC-No
#> 19 TReCASE Poisson Beta-Binomial CT1 7.07e-01 7.32e-01 7.01e-01 7.32e-01 no eQTL TReC-No;ASReC-Yes
#> 20 TReCASE Poisson Beta-Binomial CT2 0.00e+00 0.00e+00 8.47e-02 0.00e+00 cis eQTL TReC-No;ASReC-Yes
#> 21 TReCASE Poisson Beta-Binomial CT3 0.00e+00 0.00e+00 5.89e-02 0.00e+00 cis eQTL TReC-No;ASReC-Yes
#> 22 TReCASE Poisson Binomial CT1 7.07e-01 8.27e-01 4.50e-01 8.27e-01 no eQTL TReC-No;ASReC-No
#> 23 TReCASE Poisson Binomial CT2 0.00e+00 0.00e+00 2.08e-07 0.00e+00 trans eQTL TReC-No;ASReC-No
#> 24 TReCASE Poisson Binomial CT3 0.00e+00 0.00e+00 5.25e-05 0.00e+00 trans eQTL TReC-No;ASReC-NoFrom above, we see that modeling cell types through proportion and reference allele expression and correctly specifying the distribution of TReC and ASReC leads to no association or no cell type specific eQTLs.
OLS
We benchmark CSeQTL against OLS, an ordinary least squares model. This approach corresponds to Matrix eQTL (Shabalin, 2012).
Bulk eQTL
For a bulk eQTL, the model adjusts for genotype and confounders. The code below fits and tests for a bulk eQTL using the simulated dataset.
ols_out = CSeQTL_linearTest(input = sim$dat,XX = sim$XX,
RHO = sim$true_RHO,SNP = sim$true_SNP,MARG = TRUE)
names(ols_out)
#> [1] "lm_out" "out_df" "res_trim" "cooksd" "prop_trim"
# Model estimates
summary(ols_out$lm_out)
#>
#> Call:
#> lm(formula = YY_final ~ ., data = fin_XX)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.21073 -0.32724 0.01604 0.40286 1.53527
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 1.37255 0.07329 18.728 <2e-16 ***
#> X2 0.29801 0.03164 9.418 <2e-16 ***
#> X3 -0.74058 0.06288 -11.777 <2e-16 ***
#> X4 0.32274 0.03166 10.193 <2e-16 ***
#> X5 -0.33828 0.03149 -10.742 <2e-16 ***
#> SNP -0.75528 0.04456 -16.950 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.5421 on 294 degrees of freedom
#> Multiple R-squared: 0.7023, Adjusted R-squared: 0.6973
#> F-statistic: 138.7 on 5 and 294 DF, p-value: < 2.2e-16
# Effect sizes and hypothesis test
ols_out$out_df
#> CT EST SE PVAL
#> 1 Bulk -0.755285 0.04455894 0Similar to CSeQTL’s bulk testing, we have a false positive bulk eQTL.
Cell type-specific eQTLs
The code below will test for cell type-specific eQTLs with OLS.
ols_out = CSeQTL_linearTest(input = sim$dat,XX = sim$XX,
RHO = sim$true_RHO,SNP = sim$true_SNP,MARG = FALSE)
names(ols_out)
#> [1] "lm_out" "out_df" "res_trim" "cooksd" "prop_trim"
# Model estimates
summary(ols_out$lm_out)
#>
#> Call:
#> lm(formula = YY_final ~ ., data = fin_XX)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.06584 -0.29634 0.02862 0.32315 1.13060
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.09094 0.29950 0.304 0.7616
#> X2 0.30044 0.02918 10.295 <2e-16 ***
#> X3 -0.74198 0.05813 -12.764 <2e-16 ***
#> X4 0.32697 0.02922 11.191 <2e-16 ***
#> X5 -0.32246 0.02913 -11.069 <2e-16 ***
#> SNP -0.19956 0.15931 -1.253 0.2113
#> CT2 1.27072 0.32218 3.944 0.0001 ***
#> CT3 0.29099 0.42680 0.682 0.4959
#> SNP_CT2 0.28175 0.18197 1.548 0.1226
#> SNP_CT3 -0.08092 0.23251 -0.348 0.7281
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.4995 on 290 degrees of freedom
#> Multiple R-squared: 0.7508, Adjusted R-squared: 0.743
#> F-statistic: 97.07 on 9 and 290 DF, p-value: < 2.2e-16
# Effect sizes and hypothesis tests
ols_out$out_df
#> CT EST SE PVAL
#> 1 CT1 -0.19955716 0.1593063 0.21133755
#> 2 CT2 0.08219402 0.1352272 0.54378148
#> 3 CT3 -0.28047913 0.1633836 0.08710334Future Sections to create
- Trimming
- Sample data for eQTL mapping per gene and multiple SNPs
Session Info
sessionInfo()
#> R version 4.2.2 (2022-10-31 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19045)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=C
#> [2] LC_CTYPE=English_United States.utf8
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.utf8
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] CSeQTL_1.0.0 smarter_1.0.1 ggplot2_3.4.0
#>
#> loaded via a namespace (and not attached):
#> [1] BiocFileCache_2.6.0 plyr_1.8.8
#> [3] splines_4.2.2 BiocParallel_1.32.5
#> [5] TH.data_1.1-1 usethis_2.1.6
#> [7] GenomeInfoDb_1.34.9 digest_0.6.31
#> [9] htmltools_0.5.4 fansi_1.0.4
#> [11] magrittr_2.0.3 memoise_2.0.1
#> [13] remotes_2.4.2 Biostrings_2.66.0
#> [15] matrixStats_0.63.0 R.utils_2.12.2
#> [17] sandwich_3.0-2 bdsmatrix_1.3-6
#> [19] prettyunits_1.1.1 colorspace_2.1-0
#> [21] blob_1.2.4 rappdirs_0.3.3
#> [23] xfun_0.42 dplyr_1.1.2
#> [25] callr_3.7.3 crayon_1.5.2
#> [27] RCurl_1.98-1.12 jsonlite_1.8.5
#> [29] survival_3.4-0 zoo_1.8-12
#> [31] glue_1.6.2 gtable_0.3.4
#> [33] HelpersMG_5.8 zlibbioc_1.44.0
#> [35] XVector_0.38.0 DelayedArray_0.24.0
#> [37] pkgbuild_1.4.0 MatrixEQTL_2.3
#> [39] BiocGenerics_0.44.0 scales_1.2.1
#> [41] mvtnorm_1.1-3 DBI_1.1.3
#> [43] miniUI_0.1.1.1 Rcpp_1.0.10
#> [45] xtable_1.8-4 progress_1.2.2
#> [47] emdbook_1.3.12 bit_4.0.5
#> [49] stats4_4.2.2 profvis_0.3.7
#> [51] htmlwidgets_1.6.1 httr_1.4.6
#> [53] gplots_3.1.3 ellipsis_0.3.2
#> [55] farver_2.1.1 urlchecker_1.0.1
#> [57] pkgconfig_2.0.3 XML_3.99-0.14
#> [59] R.methodsS3_1.8.2 sass_0.4.5
#> [61] dbplyr_2.3.0 utf8_1.2.3
#> [63] labeling_0.4.3 tidyselect_1.2.0
#> [65] rlang_1.1.1 later_1.3.0
#> [67] AnnotationDbi_1.60.0 munsell_0.5.0
#> [69] tools_4.2.2 cachem_1.0.8
#> [71] cli_3.6.1 generics_0.1.3
#> [73] RSQLite_2.3.1 devtools_2.4.5
#> [75] evaluate_0.20 stringr_1.5.0
#> [77] fastmap_1.1.1 yaml_2.3.7
#> [79] processx_3.8.0 knitr_1.42
#> [81] bit64_4.0.5 fs_1.5.2
#> [83] caTools_1.18.2 purrr_1.0.1
#> [85] KEGGREST_1.38.0 mime_0.12
#> [87] R.oo_1.25.0 xml2_1.3.4
#> [89] biomaRt_2.54.0 compiler_4.2.2
#> [91] filelock_1.0.2 curl_5.0.1
#> [93] png_0.1-8 tibble_3.2.1
#> [95] bslib_0.4.2 stringi_1.7.12
#> [97] highr_0.10 ps_1.7.2
#> [99] GenomicFeatures_1.50.4 lattice_0.20-45
#> [101] Matrix_1.5-1 vctrs_0.6.2
#> [103] pillar_1.9.0 lifecycle_1.0.3
#> [105] jquerylib_0.1.4 data.table_1.14.8
#> [107] bitops_1.0-7 httpuv_1.6.7
#> [109] rtracklayer_1.58.0 GenomicRanges_1.50.2
#> [111] R6_2.5.1 BiocIO_1.8.0
#> [113] promises_1.2.0.1 KernSmooth_2.23-20
#> [115] IRanges_2.32.0 sessioninfo_1.2.2
#> [117] codetools_0.2-18 MASS_7.3-58.1
#> [119] gtools_3.9.4 assertthat_0.2.1
#> [121] pkgload_1.3.2 SummarizedExperiment_1.28.0
#> [123] rjson_0.2.21 withr_2.5.0
#> [125] GenomicAlignments_1.34.0 Rsamtools_2.14.0
#> [127] multcomp_1.4-20 S4Vectors_0.36.1
#> [129] GenomeInfoDbData_1.2.9 parallel_4.2.2
#> [131] hms_1.1.2 grid_4.2.2
#> [133] coda_0.19-4 rmarkdown_2.20
#> [135] MatrixGenerics_1.10.0 bbmle_1.0.25
#> [137] numDeriv_2016.8-1.1 Biobase_2.58.0
#> [139] shiny_1.7.4 restfulr_0.0.15References
Little, P., Liu, S., Zhabotynsky, V., Li, Y., Lin, D.-Y., & Sun, W.
(2023). A computational method for cell type-specific expression
quantitative trait loci mapping using bulk RNA-seq data. Nature
Communications, 14(1), 3030.
Shabalin, A. A. (2012). Matrix eQTL: Ultra fast eQTL analysis via large
matrix operations. Bioinformatics, 28(10), 1353–1358.
Sun, W. (2012). A statistical framework for eQTL mapping using RNA-seq
data. Biometrics, 68(1), 1–11.