| Title: | Anchor-Restricted Loop QTL Mapping from Hi-C |
|---|---|
| Description: | Maps chromatin loop QTLs by testing whether genotype at a loop's anchors associates with the loop's contact strength across samples. Reads .hic files via strawr, builds an exact-coordinate consensus loop set, normalizes contact strength across samples via band-wise quantile normalization in the style of BNBC (implemented on top of preprocessCore), and runs anchor-restricted SNP-by-loop regressions with MatrixEQTL. Intra-chromosomal loops only; no FDR or lead-SNP selection (left to the user). |
| Authors: | Rhie Lab at the University of Southern California [cre], Zexun Wu [aut] |
| Maintainer: | Rhie Lab at the University of Southern California <[email protected]> |
| License: | GPL-3 + file LICENSE |
| Version: | 0.99.0 |
| Built: | 2026-09-19 06:08:49 UTC |
| Source: | https://github.com/BiocStaging/loopQTL |
loopQTL objects should normally be generated by loopQTLsetup(). See its
documentation for more information on the expected format of input data.
## S4 method for signature 'loopQTL' samples(object) ## S4 method for signature 'loopQTL' loops(object) ## S4 method for signature 'loopQTL' cis(object) ## S4 method for signature 'loopQTL' vcfPath(object) ## S4 method for signature 'loopQTL' resolution(object) ## S4 method for signature 'loopQTL' phenotype(object) ## S4 method for signature 'loopQTL' results(object) ## S4 method for signature 'loopQTL' metadata(object) ## S4 method for signature 'loopQTL' show(object)## S4 method for signature 'loopQTL' samples(object) ## S4 method for signature 'loopQTL' loops(object) ## S4 method for signature 'loopQTL' cis(object) ## S4 method for signature 'loopQTL' vcfPath(object) ## S4 method for signature 'loopQTL' resolution(object) ## S4 method for signature 'loopQTL' phenotype(object) ## S4 method for signature 'loopQTL' results(object) ## S4 method for signature 'loopQTL' metadata(object) ## S4 method for signature 'loopQTL' show(object)
object |
a loopQTL object |
All slots can also be accessed through an identically named function, i.e.,
to access results(obj), use loopQTL::results(obj).
samplesList of unique sample IDs
loopsConsensus intra-chromosomal loops as a InteractionSet::GInteractions object
cisOptional numeric vector of per-sample cis-contact totals. Must
have names matching the IDs in samples, and be in the same order as
samples.
vcfPathPath to the tabix-indexed VCF
resolutionAnalysis bin size (bp)
phenotypeMatrix containing normalized numeric phenotype data for each loop (in the rows) and each sample (in the columns). Dimensions and order has not been run.
resultsQTL result table. Empty if runQTL() has not been run.
Contains a row for each loop-SNP pair and columns loopID, snpID,
chrom, pos, ref, alt, anchor, beta, se, t, p, nUsed,
and maf. If loopQTLannotate() has been run, also contains gene,
peak (optional), and Type columns.
metadataFunction-specific metadata
loopQTLsetup
obj <- loopQTLExampleData() obj samples(obj) loops(obj)obj <- loopQTLExampleData() obj samples(obj) loops(obj)
This function annotates SNP-loop QTL results with target genes, optional peak overlap, and a SNP-loop Type classification.
loopQTLannotate( obj, gtf, peak = NULL, TSSwindow = 2000, anchorWindow = 5000, peakAnnotation = "Peak", removeTSS = TRUE )loopQTLannotate( obj, gtf, peak = NULL, TSSwindow = 2000, anchorWindow = 5000, peakAnnotation = "Peak", removeTSS = TRUE )
obj |
A loopQTL object with |
gtf |
Path to a GTF file with |
peak |
Optional path to a peak bed file. If supplied, an extra |
TSSwindow |
Base-pair window flank around each gene's TSS to define promoters. Default 2000. |
anchorWindow |
Base-pair window flank added around each anchor when looking for overlapping promoters and peaks. Default 5000. |
peakAnnotation |
Character label substituted into the |
removeTSS |
Logical; default |
Adds three columns to results(obj):
geneThe name of the genes found at the target anchor. Empty if no gene is found.
peakTRUE if the SNP
overlaps any peak from the supplied peak file, FALSE otherwise. When
removeTSS = TRUE, peaks that overlap any TSS region are dropped
before this overlap check (so promoter-like peaks don't count).
TypeOne of:
"WithinTSS" - SNP overlaps any gene's TSS (± TSSwindow).
Short-circuits all other annotation when removeTSS = TRUE.
"rs-gene" - SNP not in TSS, has a target gene at the partner anchor.
"rs-<peakAnnotation>" - SNP not in TSS, no gene at partner,
but partner anchor overlaps a peak (only when peak is supplied).
e.g. "rs-Enhancer".
"rs-O" - SNP not in TSS, no gene AND no peak at the partner anchor.
"rs(<peakAnnotation>)-gene" - additionally overlaps a peak
(only when peak is supplied), has a target gene. e.g.
"rs(Enhancer)-gene".
"rs(<peakAnnotation>)-<peakAnnotation>" - SNP overlaps peak,
partner overlaps peak, no gene at partner.
e.g. "rs(Enhancer)-Enhancer".
"rs(<peakAnnotation>)-O" - overlaps a peak but no target gene
and no peak at partner.
If both promtoer and peaks are found at the target anchor, it give priority on genes.
The input obj with results(obj) augmented by gene, peak
(when applicable), and Type columns.
## Load the example object. In real analysis, the input ## should be the output of loopQTLsetup(). ##Here we just use the example object. obj <- loopQTLExampleData() ## get the path of example annotation files extdataDir <- system.file("extdata", package = "loopQTL") peakPath <- file.path(extdataDir, "example_annotation.bed") gtfPath <- file.path(extdataDir, "example.gtf.gz") ## Annotate SNP-loop pairs with gene, peak overlap, and Type obj <- loopQTLannotate( obj, gtf = gtfPath, peak = peakPath, peakAnnotation = "Enhancer", TSSwindow = 2000, removeTSS = TRUE) ##check the annotated results head(results(obj))## Load the example object. In real analysis, the input ## should be the output of loopQTLsetup(). ##Here we just use the example object. obj <- loopQTLExampleData() ## get the path of example annotation files extdataDir <- system.file("extdata", package = "loopQTL") peakPath <- file.path(extdataDir, "example_annotation.bed") gtfPath <- file.path(extdataDir, "example.gtf.gz") ## Annotate SNP-loop pairs with gene, peak overlap, and Type obj <- loopQTLannotate( obj, gtf = gtfPath, peak = peakPath, peakAnnotation = "Enhancer", TSSwindow = 2000, removeTSS = TRUE) ##check the annotated results head(results(obj))
This function runs PCA on the normalized phenotype matrix and picks a number of principal components K to use as covariates in the downstream regression. This function returns a list with the chosen K and the sample × K matrix of PCs.
loopQTLchooseK( obj, Kmethod = c("elbow", "BE"), normal = TRUE, Kmax = NULL, ... )loopQTLchooseK( obj, Kmethod = c("elbow", "BE"), normal = TRUE, Kmax = NULL, ... )
obj |
A loopQTL object with |
Kmethod |
|
normal |
Apply RINT per loop before PCA. Default |
Kmax |
The maximum number of PCs returned. Default
|
... |
Forwarded to PCAForQTL's |
A list:
KChosen number of PCs (integer).
PCssamples x K matrix of principal components.
methodWhich K-selection method was used.
obj <- loopQTLExampleData() Kres <- loopQTLchooseK(obj, Kmethod = "elbow") Kres$Kobj <- loopQTLExampleData() Kres <- loopQTLchooseK(obj, Kmethod = "elbow") Kres$K
Convenience helper for @examples blocks and vignette code that need a
fully-populated loopQTL object without re-running the full pipeline.
Loads example_loopQTL_obj.rda shipped in inst/extdata/ and rewrites
the file paths baked into the object so they resolve against the
currently installed package location.
loopQTLExampleData()loopQTLExampleData()
A loopQTL object with phenotype, results,
and metadata populated.
obj <- loopQTLExampleData() objobj <- loopQTLExampleData() obj
This function normalizes contact strength across samples so that
per-loop phenotypes are directly comparable. For each chromosome it pulls
the whole-genome matrix per sample, log-CPM adjusts, then does
band-wise quantile normalization across samples per genomic-distance band.
Normalized values at the consensus-loop bin-pairs are written to
phenotype(obj).
loopQTLnormalize( obj, BPPARAM = BiocParallel::SerialParam(), saveWindows = FALSE, windowsOut = NULL, window = 1e+06, snpWindow = 0 )loopQTLnormalize( obj, BPPARAM = BiocParallel::SerialParam(), saveWindows = FALSE, windowsOut = NULL, window = 1e+06, snpWindow = 0 )
obj |
A loopQTL object created from |
BPPARAM |
A BiocParallel::BiocParallelParam for parallelization
over chromosomes. Default |
saveWindows |
Logical. If |
windowsOut |
Required when |
window |
Minimum window extent in bp around each in-loop SNP to
create |
snpWindow |
Base-pair window flank added to each consensus anchor when selecting in-loop SNPs. Default 0. |
Please note that you need to set saveWindows to TRUE to additionally output
a per-SNP contact-window matrix as .rds file. This file is needed by the
loopQTLplotSnpWindow() and loopQTLplotLocus() downstream functions.
saveWindows does not change phenotype(obj) or affect the runQTL() step.
The loopQTL obj with phenotype(obj) filled.
## Load the example object. In real analysis, the input ## should be the output of loopQTLsetup(). ##Here we just use the example object. obj <- loopQTLExampleData() ## Example code to run loopQTLnormalize(). By default it will NOT create ## the .rds contact submatrix file that downstream functions need. obj <- loopQTLnormalize( obj, BPPARAM = BiocParallel::MulticoreParam(workers = 2), snpWindow = 10000) ## Example code to run loopQTLnormalize() and also create the `.rds` ## contact submatrix file needed by downstream visualization functions. ## First, get a temporary directory in which to save the output file. tempDirectory <- tempdir() obj <- loopQTLnormalize( obj, BPPARAM = BiocParallel::MulticoreParam(workers = 2), saveWindows = TRUE, windowsOut = file.path(tempDirectory, "loop_snp_windows.rds"), window = 2000000, snpWindow = 10000)## Load the example object. In real analysis, the input ## should be the output of loopQTLsetup(). ##Here we just use the example object. obj <- loopQTLExampleData() ## Example code to run loopQTLnormalize(). By default it will NOT create ## the .rds contact submatrix file that downstream functions need. obj <- loopQTLnormalize( obj, BPPARAM = BiocParallel::MulticoreParam(workers = 2), snpWindow = 10000) ## Example code to run loopQTLnormalize() and also create the `.rds` ## contact submatrix file needed by downstream visualization functions. ## First, get a temporary directory in which to save the output file. tempDirectory <- tempdir() obj <- loopQTLnormalize( obj, BPPARAM = BiocParallel::MulticoreParam(workers = 2), saveWindows = TRUE, windowsOut = file.path(tempDirectory, "loop_snp_windows.rds"), window = 2000000, snpWindow = 10000)
This function returns the per-sample phenotype and genotype (dosage 0/1/2) for a single SNP-loop pair. Useful for inspecting a specific hit or exporting sample-level values for reporting.
loopQTLpairData(obj, loopID, snpID, normal = TRUE)loopQTLpairData(obj, loopID, snpID, normal = TRUE)
obj |
A loopQTL object with |
loopID |
Character vector of loopIDs to extract. |
snpID |
Character vector of snpIds to extract. |
normal |
Apply RINT per loop before returning phenotype values. default |
data.frame with columns sample, loopID, snpID, genotype
(0 - Reference homozygou; 1 - heterozygous, 2 - Alternative homozygous), phenotype, cis.
## Load the example object. In real analysis, the input ## should be the output of runQTL(). ##Here we just use the example object. obj <- loopQTLExampleData() ##Check results so that you can choose your interested loopID and snpID head(results(obj)) ##Then you can specify the specifi SNP-loop pair data you want to check at samples level. pair <- loopQTLpairData( obj, loopID = "L0000187_chr1_46030000_46220000", snpID = "rs11584814") head(pair)## Load the example object. In real analysis, the input ## should be the output of runQTL(). ##Here we just use the example object. obj <- loopQTLExampleData() ##Check results so that you can choose your interested loopID and snpID head(results(obj)) ##Then you can specify the specifi SNP-loop pair data you want to check at samples level. pair <- loopQTLpairData( obj, loopID = "L0000187_chr1_46030000_46220000", snpID = "rs11584814") head(pair)
This function renders a locus-scale ggplot combining the consensus
loops, the tested SNP, and (when a GTF is supplied) gene annotations, over
a user-specified region. Please note that to use values = "cpm" or
values = "raw" to show per-loop contact strength, you need to set
saveWindows = TRUE to save the per-SNP contact-window matrices from
loopQTLnormalize().
loopQTLplotLocus( obj, snpID, chrom, start, end, gtf, pThreshold = 0.05, values = c("normalized", "cpm", "raw"), windows = NULL, strongColor = "#FF7917", weakColor = "grey90", limit = NULL, geneBiotypes = "protein_coding", axisUnit = c("Mb", "bp") )loopQTLplotLocus( obj, snpID, chrom, start, end, gtf, pThreshold = 0.05, values = c("normalized", "cpm", "raw"), windows = NULL, strongColor = "#FF7917", weakColor = "grey90", limit = NULL, geneBiotypes = "protein_coding", axisUnit = c("Mb", "bp") )
obj |
A loopQTL object with |
snpID |
Single snpID present in |
chrom |
Chromosome. e.g. "chr1". |
start |
The start genomic coordinate in bp for the x-axis. |
end |
The end genomic coordinate in bp for the x-axis. |
gtf |
Path to a GTF file. |
pThreshold |
Nominal p-value cutoff for "significant" loops. Default |
values |
Which contact-strength matrix to color arcs by:
|
windows |
Optional path/loaded list from |
strongColor |
Color for the strongest signal on the arc gradient.
Default |
weakColor |
Color for the weakest signal. Default |
limit |
Optional numeric cap on the color scale (values above are
limited to |
geneBiotypes |
Character vector of |
axisUnit |
|
A ggplot (gene / SNP / arc panels
stacked) object. Save with ggplot2::ggsave().
## Load the example object. In real analysis, the input ## should be the output of runQTL(). ## If you specify `values = "cpm"` or `values = "raw"`, the object must ## carry the contact matrix `.rds` (see `obj@metadata$windowsOut`); ## otherwise the normalized matrix stored in the object is used. ##Here we just use the example object. obj <- loopQTLExampleData() ##You also need a gtf file input to give the genomic annotations. Here we provide an example gtf. gtfPath <- system.file("extdata", "example.gtf.gz", package = "loopQTL") ## Return is a ggplot object so you can ggsave() to output. p <- loopQTLplotLocus( obj, chrom = "chr1", start = 45900000, end = 46300000, snpID = "rs11584814", values = "cpm", gtf = gtfPath)## Load the example object. In real analysis, the input ## should be the output of runQTL(). ## If you specify `values = "cpm"` or `values = "raw"`, the object must ## carry the contact matrix `.rds` (see `obj@metadata$windowsOut`); ## otherwise the normalized matrix stored in the object is used. ##Here we just use the example object. obj <- loopQTLExampleData() ##You also need a gtf file input to give the genomic annotations. Here we provide an example gtf. gtfPath <- system.file("extdata", "example.gtf.gz", package = "loopQTL") ## Return is a ggplot object so you can ggsave() to output. p <- loopQTLplotLocus( obj, chrom = "chr1", start = 45900000, end = 46300000, snpID = "rs11584814", values = "cpm", gtf = gtfPath)
This function renders the per-sample boxplot for a chosen result SNP-loop
pair. It is returned as a ggplot object that can be saved with
ggplot2::ggsave().
loopQTLplotPair(obj, loopID, snpID, normal = TRUE, jitter = TRUE)loopQTLplotPair(obj, loopID, snpID, normal = TRUE, jitter = TRUE)
obj |
A loopQTL object with |
loopID |
A single loopID to plot. |
snpID |
A single snpID to plot. |
normal |
Plot Rank-based Inverse Normal Transformed (RINT'd) phenotype
values. Defaults to TRUE. Set to |
jitter |
Add a jittered point overlay (default |
A ggplot object.
## Load the example object. In real analysis, the input ## should be the output of runQTL(). ##Here we just use the example object. obj <- loopQTLExampleData() ##Check results so that you can choose your interested loopID and snpID head(results(obj)) ##Then you can specify the specifi SNP-loop pair data you want to plot as classical QTL box plot. ##Return is a ggplot object so you can save with ggsave(). p <- loopQTLplotPair( obj, loopID = "L0000187_chr1_46030000_46220000", snpID = "rs11584814")## Load the example object. In real analysis, the input ## should be the output of runQTL(). ##Here we just use the example object. obj <- loopQTLExampleData() ##Check results so that you can choose your interested loopID and snpID head(results(obj)) ##Then you can specify the specifi SNP-loop pair data you want to plot as classical QTL box plot. ##Return is a ggplot object so you can save with ggsave(). p <- loopQTLplotPair( obj, loopID = "L0000187_chr1_46030000_46220000", snpID = "rs11584814")
This function renders a Hi-C-style heatmap centered on a chosen
SNP: for every bin-pair in the SNP's cached contact window (from
loopQTLnormalize(saveWindows = TRUE)), the color encodes the per-sample
regression coefficient (beta) of the bin-pair's contact strength on the SNP's
dosage. Set normal = TRUE and covariates = "auto" to match the model
used in runQTL().
loopQTLplotSnpWindow( obj, snpID, chrom = NULL, start = NULL, end = NULL, windows = NULL, mirror = TRUE, minSamplesPerCell = 3, normal = FALSE, covariates = NULL, Kmethod = c("elbow", "BE"), Kmax = NULL, limit = NULL, palette = "RdBu", metric = c("betaLogp", "beta"), values = c("normalized", "cpm", "raw"), axisUnit = c("Mb", "bp") )loopQTLplotSnpWindow( obj, snpID, chrom = NULL, start = NULL, end = NULL, windows = NULL, mirror = TRUE, minSamplesPerCell = 3, normal = FALSE, covariates = NULL, Kmethod = c("elbow", "BE"), Kmax = NULL, limit = NULL, palette = "RdBu", metric = c("betaLogp", "beta"), values = c("normalized", "cpm", "raw"), axisUnit = c("Mb", "bp") )
obj |
A loopQTL object with |
snpID |
A single snpID present in the saved windows file. |
chrom |
Optional chrom for verification. If supplied, must match the SNP's chrom in the windows file. |
start |
Optional Start genomic position in bp to crop the heatmap to. |
end |
Optional End genomic position in bp to crop the heatmap to. |
windows |
Either |
mirror |
Logical; mirror the upper-triangle beta values to the
lower triangle so the heatmap is symmetric. Default |
minSamplesPerCell |
Minimum usable (non-NA dose + complete covariates) samples required to fit the regression for a given cell. Cells with fewer samples are returned NA. Default 3. |
normal |
Apply Rank-based Inverse Normal Transformation (RINT) per
cell across samples before fitting. Matches the |
covariates |
Covariate handling for the regression:
|
Kmethod |
|
Kmax |
Hard cap on the auto-chosen K when |
limit |
Symmetric color-scale limit on beta. Values beyond
|
palette |
RColorBrewer palette name for the color pattern in the heatmap. Default |
metric |
What to render per cell:
|
values |
Which cell values to run the regression on:
|
axisUnit |
Unit for the graph axes - |
A ggplot object.
## Load the example object. In real analysis, the input ## should be the output of runQTL(). ##Also please note that this function requires the object to have the contact matrix `rda` file. ##Check `obj@windowsOut` to see if your object contains this file. ##Here we just use the example object. obj <- loopQTLExampleData() ##Check the results to see if you have any interested SNPs that you want to plot on. head(results(obj)) p <- loopQTLplotSnpWindow( obj, chrom = "chr1", start = 45900000, end = 46300000, snpID = "rs11584814", normal = TRUE, covariates = "auto")## Load the example object. In real analysis, the input ## should be the output of runQTL(). ##Also please note that this function requires the object to have the contact matrix `rda` file. ##Check `obj@windowsOut` to see if your object contains this file. ##Here we just use the example object. obj <- loopQTLExampleData() ##Check the results to see if you have any interested SNPs that you want to plot on. head(results(obj)) p <- loopQTLplotSnpWindow( obj, chrom = "chr1", start = 45900000, end = 46300000, snpID = "rs11584814", normal = TRUE, covariates = "auto")
validates the input files header, including .vcf and .hic files;
reads each per-sample BEDPE and builds a consensus loop set by
exact-coordinate matching. Loops found fewer than minSamples samples will be
dropped;
Automatically computes per-sample cis-contact depth from the .hic files
when the samplesheet lacks a cis column.
loopQTLsetup( samples, vcf, resolution, minSamples = 2, BPPARAM = BiocParallel::SerialParam() )loopQTLsetup( samples, vcf, resolution, minSamples = 2, BPPARAM = BiocParallel::SerialParam() )
samples |
a data.frame with columns |
vcf |
Path to a VCF. The vcf should be bgzipped with tabix index. |
resolution |
a numeric bin size recognized in the .hic file (e.g. 10000). |
minSamples |
The minimum number of samples supporting a consensus loop. Default is 2 samples. |
BPPARAM |
A BiocParallel::BiocParallelParam that parallelizes the cis-computing straw pulls. |
The samplesheet is a tab-separated file with columns Sample, hic,
loops, and optionally cis.
A loopQTL object with samples, loops, vcfPath,
resolution, cis, and metadata populated.
phenotype is an empty placeholder until loopQTLnormalize() fills it.
## Get the necessary input example files path. extdataDir <- system.file("extdata", package = "loopQTL") vcfPath <- file.path(extdataDir, "example.vcf.gz") samplesheetPath <- file.path(extdataDir, "example_samplesheet.tsv") samples <- read.table(samplesheetPath, header = TRUE) ##modify the relative file path samples$hic <- unlist( lapply(samples$hic, function(x) { file.path(extdataDir, x) })) samples$loops <- unlist( lapply(samples$loops, function(x) { file.path(extdataDir, x) })) ## Run loopQTLsetup obj <- loopQTLsetup( samples = samples, vcf = vcfPath, resolution = 10000, minSamples = 2)## Get the necessary input example files path. extdataDir <- system.file("extdata", package = "loopQTL") vcfPath <- file.path(extdataDir, "example.vcf.gz") samplesheetPath <- file.path(extdataDir, "example_samplesheet.tsv") samples <- read.table(samplesheetPath, header = TRUE) ##modify the relative file path samples$hic <- unlist( lapply(samples$hic, function(x) { file.path(extdataDir, x) })) samples$loops <- unlist( lapply(samples$loops, function(x) { file.path(extdataDir, x) })) ## Run loopQTLsetup obj <- loopQTLsetup( samples = samples, vcf = vcfPath, resolution = 10000, minSamples = 2)
This function first reads all anchor-window SNPs from the VCF, and runs
MatrixEQTL twice - once treating anchor A as the gene position, once anchor
B. The two result tables are outer-joined on SNP-loop level. Covariates are
supplied via the covariates argument - typically the phenotype PCs chosen
in the previous step. By default, rank-based inverse normal transformation
(RINT) is applied per loop before regression.
runQTL( obj, window = 0, covariates = NULL, normal = TRUE, Kmethod = c("elbow", "BE"), Kmax = NULL, pvOutputThreshold = 1 )runQTL( obj, window = 0, covariates = NULL, normal = TRUE, Kmethod = c("elbow", "BE"), Kmax = NULL, pvOutputThreshold = 1 )
obj |
A loopQTL object on which |
window |
Base-pair window flank added to each anchor bin when matching SNPs. Default 0. |
covariates |
Optional |
normal |
Apply Rank-based Inverse Normal Transformation (RINT) per loop
before regression. Default |
Kmethod |
|
Kmax |
The maximum number of PCs used. Defaults to the number of samples minus two. |
pvOutputThreshold |
MatrixEQTL p-value threshold. Default 1 (return every pair). |
The loopQTL object with results(obj) filled.
## Load the example object. In real analysis, the input ## should be the output of loopQTLnormalize(). ##Here we just use the example object. obj <- loopQTLExampleData() ##Run the function, please note that the covariates option is not specified. obj <- runQTL( obj, window = 10000) ## In real analysis we recommend running runQTL() with covariates ## specified for more accurate regression, typically the top-K PCs. ## The helper loopQTLchooseK() picks K following the PCAForQTL method ## (https://github.com/heatherjzhou/PCAForQTL). Kres <- loopQTLchooseK(obj) obj <- runQTL( obj, window = 10000, covariates = Kres$PCs)## Load the example object. In real analysis, the input ## should be the output of loopQTLnormalize(). ##Here we just use the example object. obj <- loopQTLExampleData() ##Run the function, please note that the covariates option is not specified. obj <- runQTL( obj, window = 10000) ## In real analysis we recommend running runQTL() with covariates ## specified for more accurate regression, typically the top-K PCs. ## The helper loopQTLchooseK() picks K following the PCAForQTL method ## (https://github.com/heatherjzhou/PCAForQTL). Kres <- loopQTLchooseK(obj) obj <- runQTL( obj, window = 10000, covariates = Kres$PCs)