diaQTL Vignette

diaQTL uses the Markov Chain Monte Carlo (MCMC) algorithm in R package BGLR to regress phenotypes on the expected dosage of parental origin genotypes. The example dataset is a potato half-diallel population with three founders.

Structure of the input files

Three input files are needed for QTL analysis: (1) pedigree file, (2) genotype file, (3) phenotype file. diaQTL contains several functions to prepare the input files from the output of several linkage analysis packages, including PolyOrigin (convert_polyorigin), MAPpoly (convert_mappoly), RABBIT (convert_rabbit), and onemap (convert_onemap). The tutorial dataset is a multiparental tetraploid population, for which PolyOrigin is the only option.

1) Pedigree file

The pedigree file has three columns: id, parent1, and parent2 (maternal effects are not modeled).

pedcsv <- system.file("vignette_data", "potato_ped.csv", package = "diaQTL")
ped <- read.csv(pedcsv, as.is = T)
head(ped)

##          id  parent1      parent2
## 1 W15268-1R W6511-1R VillettaRose
## 2 W15268-4R W6511-1R VillettaRose
## 3 W15268-5R W6511-1R VillettaRose
## 4 W15268-6R W6511-1R VillettaRose
## 5 W15268-7R W6511-1R VillettaRose
## 6 W15268-8R W6511-1R VillettaRose

table(apply(ped[,2:3],1,paste,collapse=" x "))

## 
## VillettaRose x W9914-1R W6511-1R x VillettaRose     W6511-1R x W9914-1R 
##                     113                     155                     147

2) Genotype file

The first 3 columns of the genotype file must be labeled marker, chrom, and cM, and the position in a reference genome (labeled bp) is optional as the fourth column (plotting features can use either cM or bp). Subsequent columns contain the genotype probabilities for each individual.

genocsv <- system.file( "vignette_data", "potato_geno.csv", package = "diaQTL" )
geno <- read.csv( genocsv, as.is = T, check.names = F )
geno[1:5,1:4]

##                marker chrom   cM     bp
## 1 solcap_snp_c2_36608     1 0.00 251132
## 2 solcap_snp_c2_36658     1 0.15 269400
## 3 solcap_snp_c1_10930     1 0.29 309342
## 4       PotVar0120130     1 0.41 353979
## 5       PotVar0120070     1 0.69 433801

Genotype probabilities are encoded as strings, following the format exported by the PolyOrigin software:

geno[1,10]

## [1] "12|13|19|23|31|34|62|63|67|69|83=>0.005|0.074|0.004|0.004|0.023|0.006|0.016|0.85|0.002|0.014|0.001"

The integers separated by | on the left side of the equal sign refer to genotype states, and the decimal numbers on the right side of the equal sign are probabilities. Only nonzero probabilities need to be included. There are 100 possible states for F1 populations, and 35 possible states for S1 populations:

library(diaQTL)
head(F1codes)

##   Code   State
## 1    1 1-1-5-5
## 2    2 1-1-5-6
## 3    3 1-1-5-7
## 4    4 1-1-5-8
## 5    5 1-1-6-6
## 6    6 1-1-6-7

head(S1codes)

##   Code   State
## 1    1 1-1-1-1
## 2    2 1-1-1-2
## 3    3 1-1-1-3
## 4    4 1-1-1-4
## 5    5 1-1-2-2
## 6    6 1-1-2-3

Each state has four integers, separated by dashes, to indicate which parental chromosomes were inherited. For F1 populations, the maternal chromosomes are labeled 1-4 and the paternal chromosomes 5-8.

3) Phenotype file

In the phenotype input file, the first column should be the individual identifier, followed by columns for different traits, and then optionally any columns with fixed effects to include in the linear model (e.g., block, environment). Only one trait, tuber shape, is provided in the example potato data set.

phenocsv <- system.file( "vignette_data", "potato_pheno.csv", package = "diaQTL" )
pheno <- read.csv( phenocsv, as.is = T )
head( pheno )

##          id tuber_shape
## 1 W15268-1R       -0.38
## 2 W15268-4R       -2.31
## 3 W15268-5R       -1.91
## 4 W15268-6R       -0.69
## 5 W15268-7R       -1.61
## 6 W15268-8R       -2.14

hist(pheno$tuber_shape,main="",xlab="Tuber shape")

To improve normality of the residuals, tuber shape in this data set is defined as log(L/W - 1), where L/W is the average length/width ratio of tubers weighing 6-10 ounces (170-285g).

Read the data

After installing and attaching the package, use read_data to read in all three files. (If there are fixed effects in the phenotype input file, they need to be specified as well; consult the reference manual.) By default, markers with the same map position in cM (using whatever numerical precision is present in the input map) are binned to reduce the computing time. The argument n.core = 2 is used for parallel execution on multiple cores.

data <- read_data(genofile = genocsv,
                  ploidy = 4,
                  pedfile = pedcsv,
                  phenofile = phenocsv, 
                  n.core = 2)

## 415 individuals with pedigree and genotype data 
## 413 individuals with pedigree, genotype and phenotype data 
## Preparing genotype data...
## Preparing for polygenic effects...

Setting parameters

The function set_params determines the burn-in and total number of iterations using the Raftery and Lewis diagnostic from R package coda, based on a 95% probability that the estimate for quantile q of the additive effects is within the interval (q-r,q+r). For the genome scan, we have found the results based on q=0.5,r=0.1 to be adequate. Because MCMC is a stochastic process, the results will not be the same each time. Results are shown for each variance component, and the user should choose values based on the slowest to converge (i.e., the largest number of iterations).

set_params( data, trait = "tuber_shape", q=0.5, r=0.1)

##          burnIn nIter
## additive     12   438
## residual      2   104

QTL Discovery

The scan1 function performs a hypothesis test for each marker bin. By default an additive model is used, which means the predictor variables are founder haplotype dosages. The test statistic is −ΔDIC, which is the DIC (Deviance Information Criterion) for the null model relative to the QTL model. Lower values of DIC indicate a better tradeoff between model complexity and goodness-of-fit. For a single hypothesis test, DIC differences of 5 or 10 are commonly recommended for model selection, but through simulation we have shown that larger differences are needed to control the genome-wide Type I error rate. The function DIC_thresh returns the −ΔDIC threshold for a half-diallel based on the number of parents, genome size, ploidy and α. The genome size for the dataset can be obtained using get_map.

get_map(data)

##     chrom Morgans
## 1       1    1.37
## 2       2    0.98
## 3       3    1.11
## 4       4    1.13
## 5       5    0.87
## 6       6    0.94
## 7       7    0.91
## 8       8    0.93
## 9       9    1.08
## 10     10    0.93
## 11     11    0.91
## 12     12    0.94
## 13 total:   12.09

alpha.05 <- DIC_thresh(genome.size=12.1,num.parents=3,
                       ploidy=4,alpha=0.05) 
alpha.1 <- DIC_thresh(genome.size=12.1,num.parents=3,
                      ploidy=4,alpha=0.1) 

ans1 <- scan1(data = data, 
              trait = "tuber_shape",
              params = list(burnIn=50,nIter=500), 
              n.core = 2)

ans1.summary <- scan1_summary(ans1, position="bp")
ans1.summary$peaks

##                   marker chrom     cM       bp        LL     deltaDIC
## 653  solcap_snp_c2_14738     1 133.09 86579271 -364.3301  -19.6959658
## 1069       PotVar0009844     2  67.97 40242967 -366.9843  -12.0168151
## 1470  solcap_snp_c1_6427     3  61.41 45124521 -366.6467  -12.3931506
## 1769 solcap_snp_c2_45936     4  15.95  3579250 -367.0495   -8.6718760
## 2257 solcap_snp_c2_23836     5   2.52   694653 -368.6811   -9.2896406
## 2848 solcap_snp_c1_15726     6   0.87   238858 -370.1426   -4.7560334
## 3499 solcap_snp_c2_19746     7  63.67 49221746 -370.9993   -4.4681446
## 3666 solcap_snp_c1_14165     8  26.70  7837521 -368.2701   -9.2643861
## 4074 solcap_snp_c2_39216     9  15.49  2792921 -370.6307   -3.5362051
## 4513       PotVar0111687    10  62.98 48721966 -291.5350 -159.9208998
## 4676 solcap_snp_c2_13392    11   0.78   597956 -373.8052   -0.1560629
## 5331  solcap_snp_c2_5463    12  93.38 60753284 -372.5445   -2.0970984

library(ggplot2)
ans1.summary$plot + geom_hline(yintercept=alpha.05,color="gold2",linetype=2) +
  geom_hline(yintercept=alpha.1,color="red",linetype=2)

The scan1_summary function returns the marker with the highest −ΔDIC score on each chromosome and a plot of the −ΔDIC profile. The peak on chromsome 10@63 cM coincides with the location of the classical Ro (round) QTL in potato, which has been identified as the gene StOFP20 (Wu et al. 2018). The 90% Bayesian credible interval (CI) can be obtained using BayesCI, based on the profile log-likelihood (LL), and this interval contains StOFP20 based on reference genome coordinates.

BayesCI(ans1,data,chrom="10",CI.prob=0.9)

##                   marker chrom    cM       bp        LL  deltaDIC
## 4509 solcap_snp_c2_45612    10 61.96 48203384 -304.3490 -132.5571
## 4510 solcap_snp_c2_49532    10 62.18 48443334 -302.4626 -136.0124
## 4511 solcap_snp_c2_25522    10 62.71 48617457 -291.1298 -158.2120
## 4512       PotVar0111667    10 62.85 48717669 -291.0995 -159.6922
## 4513       PotVar0111687    10 62.98 48721966 -291.5350 -159.9209
## 4514 solcap_snp_c2_25485    10 63.78 48737840 -294.1614 -152.9206
## 4515 solcap_snp_c2_27829    10 74.54 50782097 -308.4914 -125.1990

The small peak on chromosome 1@ 133 cM is right at the detection threshold for α = 0.05 but significant at α = 0.1.

Dominance

For a single tetraploid locus, there are four types of genetic effects. As mentioned already, additive effects are the regression coefficients for parental haplotype dosage. Digenic dominance effects are the regression coefficients for parental diplotypes, i.e., a combination of two parental haplotypes. Trigenic and quadrigenic dominance effects are also possible (and do not exist in diploid species). Throughout the diaQTL package, the argument dominance is used to specify the highest order of the effect to include in the model. Thus, dominance = 1 indicates only additive effects, dominance = 2 indicates additive and digenic effects, etc.

Having discovered a QTL, one may wish to rescan the chromosome with a digenic model to refine its position. In this case, the location of the QTL peak is unchanged, but we can see from the deltaDIC output that including digenic dominance lowered the DIC by over 15 points.

ans2 <- scan1(data = data,
              trait = "tuber_shape",
              params = list(burnIn=50,nIter=500),
              dominance = 2,
              chrom = "10",
              n.core = 2)
scan1_summary(ans2, position="bp")$peaks

##             marker chrom    cM       bp        LL  deltaDIC
## 4513 PotVar0111687    10 62.98 48721966 -274.1575 -175.1938

Markers can also be used as covariates with scan1, which is particularly useful for resolving multiple QTL on the same chromosome.

QTL Modeling

After discovery, the next step is to fit a multiple QTL model with function fitQTL. The argument qtl is a data frame with the marker names and dominance values for each QTL, and epistasis is a data frame with two markers for additive x additive epistasis. A set of progressively more complex models can be compared based on DIC. To improve accuracy, we will use more iterations based on the arguments q=0.05,r=0.025 for set_params.

qtl.10at63 <- ans1.summary$peaks$marker[10]
qtl.1at133 <- ans1.summary$peaks$marker[1]
model1 <- data.frame(marker=c(qtl.10at63,qtl.1at133),dominance=c(2,1))
model2 <- data.frame(marker=c(qtl.10at63,qtl.1at133),dominance=c(3,1))
model3 <- data.frame(marker=c(qtl.10at63,qtl.1at133),dominance=c(2,2))

set_params(data, trait = "tuber_shape", q=0.05, r=0.025, qtl = model1)

##                              burnIn nIter
## PotVar0111687 additive           14  1195
## PotVar0111687 digenic            10   982
## solcap_snp_c2_14738 additive      5   439
## residual                          3   272

params <- list(burnIn=100,nIter=5000)

fit1 <- fitQTL(data=data, trait="tuber_shape", params=params, qtl=model1)
fit2 <- fitQTL(data=data, trait="tuber_shape", params=params, qtl=model2)
fit3 <- fitQTL(data=data, trait="tuber_shape", params=params, qtl=model3)
fit4 <- fitQTL(data=data, trait="tuber_shape", params=params, qtl=model1,
               epistasis=data.frame(marker1=qtl.10at63,marker2=qtl.1at133))

fit1$deltaDIC

## [1] -194.4095

fit2$deltaDIC

## [1] -192.0648

fit3$deltaDIC

## [1] -189.5586

fit4$deltaDIC

## [1] -191.2624

Based on the DIC results, we select model1. Although inclusion of a polygenic effect does not have much impact on QTL mapping in diallel populations, it provides potentially useful information for genomic selection. The proportion of variance for each of the effects is returned in var.

fit1.poly <- fitQTL(data=data, trait="tuber_shape", params=params, qtl=model1,
                    polygenic=TRUE)
fit1.poly$deltaDIC #accept polygenic effect based on DIC

## [1] -267.0711

fit1.poly$var

##                              Mean CI.lower CI.upper
## PotVar0111687 additive       0.25     0.14     0.36
## PotVar0111687 digenic        0.09     0.03     0.18
## solcap_snp_c2_14738 additive 0.05     0.02     0.09
## polygenic                    0.29     0.19     0.40
## residual                     0.32     0.24     0.42

Plots of the additive and digenic effects for each marker are contained as elments of the list plots. For the dominance plot, digenic effects are above the diagonal, and below the diagonal is the sum of the additive and digenic effects.

fit3$plots[[qtl.10at63]]

## $additive

## 
## $digenic

Haplotype-based Selection

From the additive effects, we can select which haplotypes are desirable. In this example, since large negative values are desirable to maintain round tuber shape, the most desirable haplotype is W6511-1R.2. The function haplo_get can be used to extract the dosage of this haplotype across the population.

haplos <- haplo_get( data = data, 
                      marker = qtl.10at63)
hist(haplos[,"W6511-1R.2"],main="",xlab="Dosage")

which(haplos[,"W6511-1R.2"] > 1.8)

## W15268-53R W15268-93R W15267-10R 
##         37         59        187

The result shows there are three individuals with two copies of the W6511-1R.2 haplotype, which is possible due to “double reduction.” This occurs when a quadrivalent forms in meiosis I and sister chromatid fragments migrate to the same pole in meiosis II. The function haplo_plot can be used to visualize the pattern of recombination between parental haplotypes.

haplo_plot( data = data, 
            id = "W15268-53R", 
            chrom = 10,
            position = "cM",
            marker = qtl.10at63)

The dark blue segment indicates two copies of the W6511-1R.2 haplotype, and the dashed vertical line shows the position of the 10@63 QTL.