scan1gen examples

Generalized linear model

As a first example, let’s suppose we wanted to fit a generalized linear model for a binary phenotype, but with a probit link rather than a logit. We could use the R function glm(), as in the following. (Note that we allow just additive covariates.)

ll_glm <-
function(pr, pheno, addcovar=NULL, ...)
{
    formula <- ifelse(is.null(pr), "pheno ~ 1", "pheno ~ pr")
    if(!is.null(addcovar)) formula <- paste(formula, "+ addcovar")

    glm_out <- glm(as.formula(formula), family=binomial(link=probit))

    -glm_out$deviance/(2*log(10)) # log10 likelihood
}

Let’s now use this function in scan1gen. First we load and prepare the iron dataset included with R/qtl2.

# load R/qtl2
library(qtl2)

# read data
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))

# insert pseudomarkers into map
iron_map <- insert_pseudomarkers(iron$gmap, step=1)

# calculate genotype probabilities
iron_probs <- calc_genoprob(iron, iron_map, error_prob=0.002)

# covariates for X chr under null
Xcovar <- get_x_covar(iron)

The iron dataset includes two phenotypes; we’ll turn the first into a binary trait, 0 or 1 according to whether below or above the median.

# create binary trait
bin_pheno <- setNames(as.numeric(iron$pheno[,1] > median(iron$pheno[,1])),
                      rownames(iron$pheno))

We run scan1gen much like scan1:

out_probit <- scan1gen(iron_probs, bin_pheno, Xcovar=Xcovar, func=ll_glm, cores=0)

Compare this to the output of scan1 with a logit model. The results are very small.

out_scan1 <- scan1(iron_probs, bin_pheno, Xcovar=Xcovar, model="binary")

par(mar=c(4.1,4.1,0.6,0.6))
plot(out_probit - out_scan1, iron_map, ylim=c(-0.05,0.05),
     ylab="LOD difference (probit - logit)")

Let’s also try using glm() with a logit link, to see whether it gives the same results as scan1(). It does.

ll_glm_logit <-
function(pr, pheno, addcovar=NULL, ...)
{
    formula <- ifelse(is.null(pr), "pheno ~ 1", "pheno ~ pr")
    if(!is.null(addcovar)) formula <- paste(formula, "+ addcovar")

    glm_out <- glm(as.formula(formula), family=binomial(link=logit))

    -glm_out$deviance/(2*log(10)) # log10 likelihood
}

out_logit <- scan1gen(iron_probs, bin_pheno, Xcovar=Xcovar, func=ll_glm_logit, cores=0)

par(mar=c(4.1,4.1,0.6,0.6))
plot(out_logit - out_scan1, iron_map, ylim=c(-0.05,0.05),
     ylab="LOD difference (glm - scan1)")

Permutation test

To perform a permutation test, to establish statistical significance, use the R/qtl2 function scan1perm. We’ve added an argument scan_func so that scan1perm can be used with general scanning functions like scan1gen (and also scan1snps; see below).

To do the permutation test, we call scan1perm with scan_func=scan1gen and func=ll_glm (which gets passed along to scan1gen). Here we will seek separate autosome and X-chromosome significance thresholds.

operm <- scan1perm(iron_probs, bin_pheno, Xcovar=Xcovar, n_perm=1000,
                   perm_Xsp=TRUE, scan_func=scan1gen, func=ll_glm,
                   chr_lengths=chr_lengths(iron_map), cores=0)

Here are the 5% significance thresholds for the autosomes and X chromosome.

summary(operm)

## Autosome LOD thresholds (1000 permutations)
##      pheno1
## 0.05   3.38
##
## X chromosome LOD thresholds (28243 permutations)
##      pheno1
## 0.05   3.73

Zero-inflated Poisson regression

As a second example, let’s do a scan with a zero-inflated Poisson regression model, using the function zeroinfl in the pscl package.

Let’s start with our fitting function:

fit0infl <-
function(probs, pheno, addcovar=NULL, ...)
{
    formula <- ifelse(is.null(probs), "pheno~1", "pheno~probs")
    if(!is.null(addcovar)) formula <- paste(formula, "+ addcovar")

    pscl::zeroinfl(as.formula(formula))$loglik/log(10)
}

To illustrate the function, let’s simulate some F₂ intercross data. Our model has two QTL, with one QTL affecting the proportion of individuals with 0 phenotype, and the other affecting the average phenotype for those with non-zero phenotype.

library(qtl)

set.seed(20260613)
data(map10)
n <- 250
f2 <- sim.cross(map10, model=rbind(c(2, 45, 1, 0), c(5, 65, 1, 1)),
                n.ind=n, type="f2")

f2$pheno[,1] <- rbinom(nind(f2), 1, p=c(0.8, 0.6, 0.4)[f2$qtlgeno[,1]])
f2$pheno[,2] <- f2$pheno[,1]
nonzero <- (f2$pheno[,1]>0)
f2$pheno[nonzero,2] <- rpois(sum(nonzero), c(6,4,2)[f2$qtlgeno[nonzero,2]])

f2$pheno$sex <- factor(rep("male", n))
colnames(f2$pheno) <- c("zero", "count", "sex")

f2 <- convert2cross2(f2)

We can run the genome scan as follows:

f2_map <- insert_pseudomarkers(f2$gmap, step=1)
f2_probs <- calc_genoprob(f2, f2_map, err=0)

out_0infl <- scan1gen(f2_probs, f2$pheno[,"count"], func=fit0infl, cores=0)

Here’s a plot of the results:

par(mar=c(5.1,4.1,0.6,0.6))
plot(out_0infl, f2_map)

We can perform a permutation test as we did with glm:

operm_0infl <- scan1perm(f2_probs, f2$pheno[,"count"], scan_func=scan1gen,
                         func=fit0infl, cores=0, n_perm=1000)

Here is the 5% significance threshold:

summary(operm_0infl)

## LOD thresholds (1000 permutations)
##      pheno1
## 0.05   4.76

regress package

As a third example, let’s try using the regress package to fit a polygenic model. While scan1 in R/qtl2 estimates the polygenic variance under the null hypothesis of no QTL, and then takes it as fixed when scanning the genome, here we will re-estimate the residual polygenic variance at each QTL location.

Here is our fitting function:

fit_regress <-
function(probs, pheno, addcovar=NULL, kinship, ...)
{
    formula <- ifelse(is.null(probs), "pheno~1", "pheno~probs")
    K <- model.matrix(~rep(1, length(pheno)))

    if(!is.null(addcovar)) formula <- paste(formula, "+ addcovar")

    regress::regress(as.formula(formula), ~kinship, kernel=K)$llik/log(10)
}

We’ll use the iron dataset again. We need to calculate the kinship matrix. We can try either the overall matrix, or use “loco” (leave one chromosome out).

k <- calc_kinship(iron_probs)
k_loco <- calc_kinship(iron_probs, "loco")

We’ll now do the genome scan four ways: with scan1 or regress, and with the overall kinship or the “loco” kinship.

out_scan1 <- scan1(iron_probs, iron$pheno, Xcovar=Xcovar, kinship=k,
                   cores=0)
out_regress <- scan1gen(iron_probs, iron$pheno, Xcovar=Xcovar,
                        func=fit_regress, kinship=k, cores=0)
out_scan1_loco <- scan1(iron_probs, iron$pheno, Xcovar=Xcovar,
                        kinship=k_loco, cores=0)
out_regress_loco <- scan1gen(iron_probs, iron$pheno, Xcovar=Xcovar,
                             func=fit_regress, kinship=k_loco, cores=0)

Here’s a plot of the four scans for the liver phenotype: The scan1 results are in blue and the regress results are in pink. Solid curves are with an overall kinship matrix, and dashed curves are with “loco”. The results using regress are hardly different; not using loco seriously dampens QTL evidence.

par(mar=c(5.1, 4.1, 0.6, 0.6))
plot(out_scan1, iron_map, ylim=c(0, 7))
plot(out_scan1_loco, iron_map, add=TRUE, lty=2)
plot(out_regress, iron_map, add=TRUE, col="violetred")
plot(out_regress_loco, iron_map, add=TRUE, lty=2, col="violetred")

Here are the results for the spleen phenotype, which show similar differences.

par(mar=c(5.1, 4.1, 0.6, 0.6))
plot(out_scan1, iron_map, lodcolumn=2, ylim=c(0, 13))
plot(out_scan1_loco, iron_map, lodcolumn=2, add=TRUE, lty=2)
plot(out_regress, iron_map, lodcolumn=2, add=TRUE, col="violetred")
plot(out_regress_loco, iron_map, lodcolumn=2, add=TRUE, lty=2, col="violetred")

It’s easier to see the differences between the scan1 and regress if we subtract the LOD scores. Here are the results with the liver phenotype.

par(mar=c(5.1, 4.1, 0.6, 0.6))
plot(out_regress - out_scan1, iron_map, ylim=c(-0.25, 0.25))
plot(out_regress_loco - out_scan1_loco, iron_map, add=TRUE, lty=2)

And here are the LOD differences for the spleen phenotype.

par(mar=c(5.1, 4.1, 0.6, 0.6))
plot(out_regress - out_scan1, iron_map, lodcolumn=2, ylim=c(-0.25, 0.25))
plot(out_regress_loco - out_scan1_loco, iron_map, lodcolumn=2, add=TRUE, lty=2)

Finally, here’s how to do a permutation test. We’ll just use the “loco” method. To save time, we’ll just do 100 permutations.

operm_regress_loco <- scan1perm(iron_probs, iron$pheno, Xcovar=Xcovar,
                                scan_func=scan1gen, func=fit_regress,
                                kinship=k_loco, n_perm=100, cores=0)

Here’s the 5% significance threshold:

summary(operm_regress_loco)

## LOD thresholds (100 permutations)
##      liver spleen
## 0.05  3.39   3.34

sommer package

As our final example of scangen`, let’s do basically the same thing as we did with the regress package, but using the sommer package.

Here’s the fitting function:

fit_sommer <- function(probs, pheno, addcovar=NULL, kinship, ...)
{
    df <- as.data.frame(pheno)
    colnames(df) <- "pheno"

    if(is.null(probs)) {
        formula <- "pheno ~ 1"
    } else {
        df <- cbind(df, probs)
        formula <- paste("pheno ~", paste(colnames(probs), collapse="+"))
    }

    if(!is.null(addcovar)) {
        formula <- paste(formula, "+", paste(colnames(addcovar), collapse="+"))
        df <- cbind(df, addcovar)
    }
    df <- cbind(df, id=rownames(pheno))

    max(sommer::mmes(as.formula(formula), random=~sommer::vsm(sommer::ism(id), Gu=kinship), data=df,
                     dateWarning=FALSE, verbose=FALSE)$llik)/log(10)
}

Let’s try using the grav2 dataset included with R/qtl2. This is for a set of Arabidopsis recombinant inbred lines, with the outcome being a gravitropism phenotype.

grav2 <- read_cross2(system.file("extdata", "grav2.zip", package="qtl2"))
grav2_map <- insert_pseudomarkers(grav2$gmap, step=1)
grav2_probs <- calc_genoprob(grav2, grav2_map, cores=0)
grav2_k <- calc_kinship(grav2_probs, "loco")

We’ll first do genome scans for all phenotypes (“tip angle” at different times), using scan1.

out_scan1 <- scan1(grav2_probs, grav2$pheno, kinship=grav2_k, cores=0)

Let’s pick out just the phenotype with maximum LOD score, and use that with the sommer package:

mxphe <- which.max(apply(out_scan1, 2, max))

Now the genome scan with our fitting function that uses sommer:

out_sommer <- scan1gen(grav2_probs, grav2$pheno[, mxphe], kinship=grav2_k,
                       func=fit_sommer, cores=0)

Here are plots of the two:

plot(out_scan1, grav2_map, lodcolumn=mxphe)
plot(out_sommer, grav2_map, add=TRUE, col="violetred")

Here is a plot of the differences:

plot(out_sommer - out_scan1[,mxphe], grav2_map, ylim=c(-0.2, 0.2))

Here is how we would do a permutation test. To save time, we’ll just do 100 permutations.

operm_sommer <- scan1perm(grav2_probs, grav2$pheno[,mxphe], kinship=grav2_k,
                          scan_func=scan1gen, func=fit_sommer,
                          n_perm=100, cores=0)

And here is the 5% significance threshold:

summary(operm_sommer)

## LOD thresholds (100 permutations)
##      pheno1
## 0.05   3.08

Treatment of multiple phenotypes

Just like scan1, the scan1gen function can handle multiple phenotypes. But we assume the user-supplied fitting function handles just a single phenotype, and we then call it repeatedly for each of the phenotype columns.

You could, instead, handle the multiple phenotypes within the fitting function. In some cases, this can be done with considerable savings in computation time. The fitting function should then return a vector of LOD scores, one for each phenotype column, and when you call scan1gen, you should use vectorize_func=FALSE.

Permutation tests with scan1snps

As a final example, we’ll show how to use scan1perm with the scan1snps function: to do permutations when performing a SNP-by-SNP association scan in a multi-parent population.

We’ll consider the Diversity Outcross data from Gatti et al. (2014), available at the qtl2data repository.

First we load the data:

file <- "https://raw.githubusercontent.com/rqtl/qtl2data/main/DO_Gatti2014/do.zip"
do <- read_cross2(file)

Let’s calculate genotype probabilities and the “loco” kinship.

do_gmap <- insert_pseudomarkers(do$gmap, step=1)
do_pmap <- interp_map(do_gmap, do$gmap, do$pmap)
do_probs <- calc_genoprob(do, do_gmap, error_prob=0.002, cores=0)
do_aprobs <- genoprob_to_alleleprob(do_probs, cores=0)
do_k <- calc_kinship(do_aprobs, "loco", cores=0)

We’ll consider just the WBC phenotype, and we’ll use a square-root transformation.

do_pheno <- sqrt(do$pheno[,1])

To perform a SNP-based genome scan, we need access to a database of SNPs in the founder lines, available at figshare. We then use create_variant_query_func to create a function that queries the database.

query_variants <- create_variant_query_func("~/Data/CCdb/cc_variants.sqlite")

To perform a SNP association scan, we use scan1snps():

out_snps <- scan1snps(do_aprobs, map=do_pmap, pheno=do_pheno, kinship=do_k,
                      query_func=query_variants, cores=0)

And to do a permutation test, we use scan1perm with scan_func=scan1snps. We’ll just do 100 permutations, to save time.

operm_snps <- scan1perm(do_aprobs, map=do_pmap, pheno=do_pheno, kinship=do_k,
                        scan_func=scan1snps, query_func=query_variants,
                        cores=0, n_perm=100)

Here’s the 5% significance threshold:

summary(operm_snps)

## LOD thresholds (100 permutations)
##      pheno1
## 0.05   5.15

Karl Broman