As detailed in this page and in this video, the following linkage map was built based in a F2 population, derived from an interspecific cross between Mimulus guttatus and Mimulus nasutus.

Interval Mapping (IM) for the phenotype corolla width

Based on the map above, an Interval Mapping (IM) approach was performed to find QTLs, as follows.

Exportation of the list with the linkage groups (markers and positions) from package OneMap

write_map(maps_list, "mimulus_onemap.map")

Installation and loading the R/qtl package (more information about it in this page)

install.packages("qtl")

library("qtl")

Directory where we can find the raw_file with genotypic and phenotypic information

raw_file <- paste(system.file("extdata", package = "onemap"),
                  "m_feb06.raw", sep = "/")

Data reading

f2_qtlmimulus <- read.cross("mm", file = raw_file, mapfile = "mimulus_onemap.map")

##  --Read the following data:
##  Type of cross:          f2 
##  Number of individuals:  287 
##  Number of markers:      418 
##  Number of phenotypes:   16 
##  --Cross type: f2

f2_qtlmimulus

##   This is an object of class "cross".
##   It is too complex to print, so we provide just this summary.
##     F2 intercross
## 
##     No. individuals:    287 
## 
##     No. phenotypes:     16 
##     Percent phenotyped: 96.2 93.4 96.2 93 95.8 87.8 96.2 95.8 96.2 96.2 
##                         90.2 96.2 95.8 96.2 95.8 96.2 
## 
##     No. chromosomes:    14 
##         Autosomes:      1 2 3 4 5 6 7 8 9 10 11 12 13 14 
## 
##     Total markers:      390 
##     No. markers:        20 22 26 35 21 46 32 26 25 21 35 24 19 38 
##     Percent genotyped:  89.6 
##     Genotypes (%):      AA:16.4  AB:24.8  BB:22.4  not BB:19.3  not AA:17.1

Marker information plots

plotMap(f2_qtlmimulus,main=substitute(paste("Genetic linkage map - F2 ")(italic('Mimulus guttatus x Mimulus nasutus'))))

plotMissing(f2_qtlmimulus, main="Missing genotypes")

Corolla width (mm)

#plot.pheno<-plot.pheno(f2_qtlmimulus, pheno.col = 16)
plot(plot.pheno, col="yellow", xlab="Corola width (mm)", main="ww")

Jitter map

mimulus.jm_im <- jittermap(f2_qtlmimulus, amount=1e-6)
summary(mimulus.jm_im)

##     F2 intercross
## 
##     No. individuals:    287 
## 
##     No. phenotypes:     16 
##     Percent phenotyped: 96.2 93.4 96.2 93 95.8 87.8 96.2 95.8 96.2 96.2 
##                         90.2 96.2 95.8 96.2 95.8 96.2 
## 
##     No. chromosomes:    14 
##         Autosomes:      1 2 3 4 5 6 7 8 9 10 11 12 13 14 
## 
##     Total markers:      390 
##     No. markers:        20 22 26 35 21 46 32 26 25 21 35 24 19 38 
##     Percent genotyped:  89.6 
##     Genotypes (%):      AA:16.4  AB:24.8  BB:22.4  not BB:19.3  not AA:17.1

Conditional genotype probability in each genome position

mimulusprob <- calc.genoprob(mimulus.jm_im, step = 1, off.end=0, error.prob=0.001, map.function="kosambi", stepwidth="fixed")

Genotype simulation given observed marker data

mimulus1 <- sim.geno(mimulusprob, n.draws=16, step=1, off.end=0, error.prob=0.001, map.function="kosambi")

Genome scan with a single QTL model
- EM algorithm method

out_emmimulus1 <- scanone(mimulus1, pheno.col=16, method = "em")

## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 11 individuals with missing phenotypes.

head(out_emmimulus1)

##         chr pos       lod
## AAT225    1   0 0.5475985
## c1.loc1   1   1 0.5106363
## c1.loc2   1   2 0.4742116
## c1.loc3   1   3 0.4405052
## c1.loc4   1   4 0.4117262
## c1.loc5   1   5 0.3899564

summary(out_emmimulus1)

##            chr   pos   lod
## c1.loc55     1  55.0 1.502
## c2.loc49     2  49.0 2.824
## c3.loc94     3  94.0 5.156
## MgSTS455     4 209.0 1.589
## c5.loc12     5  12.0 1.694
## MgSTS542A    6  33.7 2.341
## AAT296       7 164.1 1.983
## BB176        8  89.5 1.714
## BC388        9 168.4 0.433
## MgSTS622    10 132.6 1.420
## c11.loc162  11 162.0 2.601
## c12.loc60   12  60.0 2.640
## CC330       13 126.6 0.499
## c14.loc113  14 113.0 9.034

plot(out_emmimulus1, col ="blue", ylab="LOD score", xlab="Linkage Group")

Haley-Knott regression method

out_hkmimulus1 <- scanone(mimulus1, pheno.col=16, method = "hk")

## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 11 individuals with missing phenotypes.

summary(out_hkmimulus1)

##            chr   pos   lod
## c1.loc56     1  56.0 1.605
## c2.loc49     2  49.0 2.810
## c3.loc94     3  94.0 5.348
## MgSTS455     4 209.0 1.525
## c5.loc12     5  12.0 1.735
## MgSTS542A    6  33.7 2.311
## AAT296       7 164.1 1.993
## BB176        8  89.5 1.707
## c9.loc159    9 159.0 0.366
## MgSTS622    10 132.6 1.406
## c11.loc161  11 161.0 2.596
## c12.loc60   12  60.0 2.550
## CC330       13 126.6 0.570
## c14.loc113  14 113.0 9.084

plot(out_hkmimulus1, col = "red", ylab="LOD score", xlab="Linkage Group")

Plot: Interval Mapping methods comparison: EM algorithm and Haley-Knott

plot(out_emmimulus1, out_hkmimulus1, col = c("blue", "red"), ylab="LOD score", xlab="Linkage Group")

Permutation test for significant LOD threshold definition (using HK regression)

operm.hk_im <- scanone(mimulus1,  pheno.col=16, method="hk", n.perm=1000)

## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 11 individuals with missing phenotypes.

## Doing permutation in batch mode ...

threshold_hk_im <- summary(operm.hk_im)
summary(operm.hk_im,alpha=0.05)

## LOD thresholds (1000 permutations)
##     lod
## 5% 3.78

Permutation distribution plot: significant threshold in percentile 95

plot.hk_im<-plot(operm.hk_im, col="blue")
abline(plot.hk_im, v=threshold_hk_im[1,1], col="red", lty=3, lwd=3)

Plot: Interval Mapping (EM algorithm and Haley-Knott methods)

plot.hk.em<-plot(out_emmimulus1, out_hkmimulus1, ylab="LOD score", xlab="Linkage group", lty=1, col=c("blue","red"), ylim=c(0,9))
abline(plot.hk.em, h= threshold_hk_im[1,1], col="red", lty=3)

QTL detection

summary_hk_im <- summary(out_hkmimulus1, perms=operm.hk_im, alpha = 0.05)
summary_hk_im

##            chr pos  lod
## c3.loc94     3  94 5.35
## c14.loc113  14 113 9.08

QTL observation on the linkage map

qtl_im <- makeqtl(mimulusprob, chr=summary_hk_im$chr, pos=summary_hk_im$pos, what="prob")
plot(qtl_im)

QTL additive and dominance effects detection

\(a=\displaystyle\frac{(\mu_{BB} - \mu_{AA})}{2}\)

\(d\)=\(\mu_{AB} - \displaystyle\frac{(\mu_{AA} + \mu_{BB})}{2}\)

Linkage Groups 1 to 7

chr1_7_im<-effectscan(mimulus1, chr=c(1,2,3,4,5,6,7), pheno.col=16, get.se=FALSE, draw=TRUE,
gap=25, mtick=c("line","triangle"),add.legend=TRUE, alternate.chrid=FALSE, ylab="Effect", xlab="Linkage group")

summary(chr1_7_im)

##           chr pos     a       d
## c1.loc78    1  78 0.823  0.3718
## c2.loc46    2  46 0.969 -0.2659
## c3.loc103   3 103 1.247 -1.0018
## c4.loc6     4   6 0.496  0.1741
## c5.loc90    5  90 0.622  0.4467
## c6.loc5     6   5 0.559 -0.5807
## c7.loc69    7  69 0.557 -0.0579

Linkage Groups 8 to 14

chr8_14_im<-effectscan(mimulus1, chr=c(8,9,10,11,12,13,14), pheno.col=16, get.se=FALSE, draw=TRUE,
gap=25, mtick=c("line","triangle"),add.legend=TRUE, alternate.chrid=FALSE, ylab="Effect", xlab="Linkage group")

summary(chr8_14_im)

##            chr pos     a        d
## c8.loc93     8  93 0.690 -0.26763
## c9.loc39     9  39 0.390 -0.00665
## MgSTS622    10 133 0.652  0.44197
## MgSTS494    11 196 2.202 -1.69183
## c12.loc114  12 114 1.884  1.61053
## c13.loc78   13  78 0.339  0.11066
## c14.loc113  14 113 1.858  0.14005

Fitting the multiple-QTL model

out_im <- fitqtl(mimulus1, pheno.col=16, qtl=qtl_im, method="hk", get.ests=TRUE)

## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 11 individuals with missing phenotypes.

summary(out_im)

## 
##      fitqtl summary
## 
## Method: Haley-Knott regression 
## Model:  normal phenotype
## Number of observations : 276 
## 
## Full model result
## ----------------------------------  
## Model formula: y ~ Q1 + Q2 
## 
##        df        SS         MS      LOD     %var Pvalue(Chi2)    Pvalue(F)
## Model   4  561.2411 140.310267 12.25079 18.48713 1.639533e-11 2.437561e-11
## Error 271 2474.6059   9.131387                                            
## Total 275 3035.8469                                                       
## 
## 
## Drop one QTL at a time ANOVA table: 
## ----------------------------------  
##          df Type III SS   LOD  %var F value Pvalue(Chi2) Pvalue(F)    
## 3@94.0    2       134.3 3.167 4.423   7.352        0.001  0.000777 ***
## 14@113.0  2       302.1 6.903 9.951  16.542        0.000  1.67e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Estimated effects:
## -----------------
##               est      SE      t
## Intercept 15.7787  0.1870 84.381
## 3@94.0a    1.1895  0.3111  3.824
## 3@94.0d   -0.6199  0.5119 -1.211
## 14@113.0a  1.6458  0.2863  5.748
## 14@113.0d  0.1652  0.3754  0.440

Table 1: Linkage Group, position, LOD score, QTLs effects (additive and dominant) and determination coefficient (%) for the Interval Mapping (IM) for the phenotype corolla width, using the Haley-Knott regression.

table1<-data.frame("Linkage Group"=c(summary_hk_im$chr[1],summary_hk_im$chr[2]),
                     "Position"=c(summary_hk_im$pos[1],summary_hk_im$pos[2]),
                     "LOD"=c(round(summary_hk_im$lod[1],4),round(summary_hk_im$lod[2],4)),
                   "Additive effect"=c(round(summary(chr1_7_im)$a[3],4),round(summary(chr8_14_im)$a[7],4)),
                   "Dominance Effect"=c(round(summary(chr1_7_im)$d[3],4),round(summary(chr8_14_im)$d[7],4)),
                                       "R2"=c(round(summary(out_im)$result.drop[1,4],4),round(summary(out_im)$result.drop[2,4],4)))
knitr::kable(table1)

Linkage.Group	Position	LOD	Additive.effect	Dominance.Effect	R2
3	94	5.3476	1.2473	-1.0018	4.4229
14	113	9.0840	1.8578	0.1400	9.9509

Conclusion (Interval Mapping-IM): using both algorithm EM and Haley-Knott methods and considering the LOD threshold obtained through permutations, we observed two possible QTLs for the phenotype corolla width: one in Linkage Group 3 and other in Linkage Group 14, presenting mainly additive effects. As we can see by the determination coefficient, the percentage of explained phenotypic variance for the trait by the QTLs in the linkage groups 3 and 14 was 4.37 and 9.97, respectively.

Composite Interval Mapping (CIM) for the phenotype corolla width

Loading the R/qtl package (more information about it in this page)

library("qtl")

Directory where we can find the raw_file with genotypic and phenotypic information

raw_file <- paste(system.file("extdata", package = "onemap"),
                  "m_feb06.raw", sep = "/")

Data reading

f2_qtlmimulus <- read.cross("mm", file = raw_file, mapfile = "mimulus_onemap.map")

##  --Read the following data:
##  Type of cross:          f2 
##  Number of individuals:  287 
##  Number of markers:      418 
##  Number of phenotypes:   16 
##  --Cross type: f2

summary(f2_qtlmimulus)

##     F2 intercross
## 
##     No. individuals:    287 
## 
##     No. phenotypes:     16 
##     Percent phenotyped: 96.2 93.4 96.2 93 95.8 87.8 96.2 95.8 96.2 96.2 
##                         90.2 96.2 95.8 96.2 95.8 96.2 
## 
##     No. chromosomes:    14 
##         Autosomes:      1 2 3 4 5 6 7 8 9 10 11 12 13 14 
## 
##     Total markers:      390 
##     No. markers:        20 22 26 35 21 46 32 26 25 21 35 24 19 38 
##     Percent genotyped:  89.6 
##     Genotypes (%):      AA:16.4  AB:24.8  BB:22.4  not BB:19.3  not AA:17.1

Jitter map

mimulus.jm_cim <- jittermap(f2_qtlmimulus, amount=1e-6)
summary(mimulus.jm_cim)

##     F2 intercross
## 
##     No. individuals:    287 
## 
##     No. phenotypes:     16 
##     Percent phenotyped: 96.2 93.4 96.2 93 95.8 87.8 96.2 95.8 96.2 96.2 
##                         90.2 96.2 95.8 96.2 95.8 96.2 
## 
##     No. chromosomes:    14 
##         Autosomes:      1 2 3 4 5 6 7 8 9 10 11 12 13 14 
## 
##     Total markers:      390 
##     No. markers:        20 22 26 35 21 46 32 26 25 21 35 24 19 38 
##     Percent genotyped:  89.6 
##     Genotypes (%):      AA:16.4  AB:24.8  BB:22.4  not BB:19.3  not AA:17.1

Conditional genotype probability in each genome position

mimulus.prob <- calc.genoprob(mimulus.jm_cim, step=1, error.prob=0.001, map.function="kosambi", stepwidth="fixed")

Genotype simulation given observed marker data

mimulus.qtl_cim <- sim.geno(mimulus.prob, n.draws=16, step=1, off.end=0, error.prob=0.001,
map.function="kosambi", stepwidth="fixed")

Covariates selection: stepwise procedure
- Markers in the genome with some influence on the QTL detection
- Function: control the variation out of the interval under analysis: better precision
- Covariates Selection
  - few covariates: no reduction of residual variance
  - many covariates: superparametrization, reducing power analysis
  - stepwise procedure: avoid marker selection near the putative QTL

stepwiseqtl(mimulus.qtl_cim, pheno.col=16, method="hk",  model="normal", refine.locations=TRUE, verbose=TRUE)

##  -Initial scan
## initial lod:  9.084022 
## ** new best ** (pLOD increased by 5.564)
##     no.qtl =  1   pLOD = 5.564022   formula: y ~ Q1 
##  -Step 1 
##  ---Scanning for additive qtl
##         plod = 5.568528 
##  ---Scanning for QTL interacting with Q1
##         plod = 5.416851 
##  ---Refining positions
##  ---  Moved a bit
##     no.qtl =  2   pLOD = 5.598872   formula: y ~ Q1 + Q2 
## ** new best ** (pLOD increased by 0.0348)
##  -Step 2 
##  ---Scanning for additive qtl
##         plod = 4.909327 
##  ---Scanning for QTL interacting with Q1
##         plod = 4.000398 
##  ---Scanning for QTL interacting with Q2
##         plod = 4.475696 
##  ---Look for additional interactions
##         plod = 5.357846 
##  ---Refining positions
##  ---  Moved a bit
##     no.qtl =  2   pLOD = 5.645633   formula: y ~ Q1 + Q2 + Q1:Q2 
## ** new best ** (pLOD increased by 0.0468)
##  -Step 3 
##  ---Scanning for additive qtl
##         plod = 4.931868 
##  ---Scanning for QTL interacting with Q1
##         plod = 2.369091 
##  ---Scanning for QTL interacting with Q2
##         plod = 2.189114 
##  ---Look for additional interactions
##  ---Refining positions
##     no.qtl =  3   pLOD = 4.931868   formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 
##  -Step 4 
##  ---Scanning for additive qtl
##         plod = 3.617729 
##  ---Scanning for QTL interacting with Q1
##         plod = 1.470267 
##  ---Scanning for QTL interacting with Q2
##         plod = 1.31221 
##  ---Scanning for QTL interacting with Q3
##         plod = 4.167961 
##  ---Look for additional interactions
##         plod = 1.920511 
##  ---Refining positions
##  ---  Moved a bit
##     no.qtl =  4   pLOD = 4.939461   formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4 + Q3:Q4 
##  -Step 5 
##  ---Scanning for additive qtl
##         plod = 3.596101 
##  ---Scanning for QTL interacting with Q1
##         plod = 2.417985 
##  ---Scanning for QTL interacting with Q2
##         plod = 1.534324 
##  ---Scanning for QTL interacting with Q3
##         plod = 1.042154 
##  ---Scanning for QTL interacting with Q4
##         plod = 5.80342 
##  ---Look for additional interactions
##         plod = 0.304798 
##  ---Refining positions
##  ---  Moved a bit
##     no.qtl =  5   pLOD = 6.469407   formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4 + Q3:Q4 + Q5 + Q4:Q5 
## ** new best ** (pLOD increased by 0.8238)
##  -Step 6 
##  ---Scanning for additive qtl
##         plod = 5.379587 
##  ---Scanning for QTL interacting with Q1
##         plod = 2.442267 
##  ---Scanning for QTL interacting with Q2
##         plod = 3.041442 
##  ---Scanning for QTL interacting with Q3
##         plod = 3.055628 
##  ---Scanning for QTL interacting with Q4
##         plod = 3.744503 
##  ---Scanning for QTL interacting with Q5
##         plod = 6.283052 
##  ---Look for additional interactions
##         plod = 2.446438 
##  ---Refining positions
##  ---  Moved a bit
##     no.qtl =  6   pLOD = 6.600555   formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4 + Q3:Q4 + Q5 + Q4:Q5 + Q6 + Q5:Q6 
## ** new best ** (pLOD increased by 0.1311)
##  -Step 7 
##  ---Scanning for additive qtl
##         plod = 5.277868 
##  ---Scanning for QTL interacting with Q1
##         plod = 2.768624 
##  ---Scanning for QTL interacting with Q2
##         plod = 3.550894 
##  ---Scanning for QTL interacting with Q3
##         plod = 2.218762 
##  ---Scanning for QTL interacting with Q4
##         plod = 3.716098 
##  ---Scanning for QTL interacting with Q5
##         plod = 2.645031 
##  ---Scanning for QTL interacting with Q6
##         plod = 4.218946 
##  ---Look for additional interactions
##         plod = 3.137856 
##  ---Refining positions
##  ---  Moved a bit
##     no.qtl =  7   pLOD = 5.919536   formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4 + Q3:Q4 + Q5 + Q4:Q5 + Q6 + Q5:Q6 + Q7 
##  -Step 8 
##  ---Scanning for additive qtl
##         plod = 4.770994 
##  ---Scanning for QTL interacting with Q1
##         plod = 2.732202 
##  ---Scanning for QTL interacting with Q2
##         plod = 3.05775 
##  ---Scanning for QTL interacting with Q3
##         plod = 1.505382 
##  ---Scanning for QTL interacting with Q4
##         plod = 1.847855 
##  ---Scanning for QTL interacting with Q5
##         plod = 1.929922 
##  ---Scanning for QTL interacting with Q6
##         plod = 3.72334 
##  ---Scanning for QTL interacting with Q7
##         plod = 2.935511 
##  ---Look for additional interactions
##         plod = 2.939105 
##  ---Refining positions
##  ---  Moved a bit
##     no.qtl =  8   pLOD = 4.792177   formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4 + Q3:Q4 + Q5 + Q4:Q5 + Q6 + Q5:Q6 + Q7 + Q8 
##  -Step 9 
##  ---Scanning for additive qtl
##         plod = 3.675926 
##  ---Scanning for QTL interacting with Q1
##         plod = 1.247319 
##  ---Scanning for QTL interacting with Q2
##         plod = 2.573498 
##  ---Scanning for QTL interacting with Q3
##         plod = 0.8255097 
##  ---Scanning for QTL interacting with Q4
##         plod = 0.8421086 
##  ---Scanning for QTL interacting with Q5
##         plod = 0.8188575 
##  ---Scanning for QTL interacting with Q6
##         plod = 1.850082 
##  ---Scanning for QTL interacting with Q7
##         plod = 2.103274 
##  ---Scanning for QTL interacting with Q8
##         plod = 3.027938 
##  ---Look for additional interactions
##         plod = 1.87568 
##  ---Refining positions
##  ---  Moved a bit
##     no.qtl =  9   pLOD = 3.69499   formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4 + Q3:Q4 + Q5 + Q4:Q5 + Q6 + Q5:Q6 + Q7 + Q8 + Q9 
##  -Step 10 
##  ---Scanning for additive qtl
##         plod = 2.366547 
##  ---Scanning for QTL interacting with Q1
##         plod = 0.2553147 
##  ---Scanning for QTL interacting with Q2
##         plod = 1.5149 
##  ---Scanning for QTL interacting with Q3
##         plod = -0.2092309 
##  ---Scanning for QTL interacting with Q4
##         plod = 0.08071932 
##  ---Scanning for QTL interacting with Q5
##         plod = -0.0235896 
##  ---Scanning for QTL interacting with Q6
##         plod = 0.4215054 
##  ---Scanning for QTL interacting with Q7
##         plod = 1.578466 
##  ---Scanning for QTL interacting with Q8
##         plod = 2.185365 
##  ---Scanning for QTL interacting with Q9
##         plod = 4.138699 
##  ---Look for additional interactions
##         plod = 2.673088 
##  ---Refining positions
##  ---  Moved a bit
##     no.qtl =  10   pLOD = 4.357977   formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4 + Q3:Q4 + Q5 + Q4:Q5 + Q6 + Q5:Q6 + Q7 + Q8 + Q9 + Q10 + Q9:Q10 
##  -Starting backward deletion
##  ---Dropping Q8 
##     no.qtl =  9   pLOD = 5.653258   formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q7 + Q8 + Q9 + Q1:Q2 + Q3:Q4 + Q4:Q5 + Q5:Q6 + Q8:Q9 
##  ---Refining positions
##  ---  Moved a bit
##  ---Dropping Q7 
##     no.qtl =  8   pLOD = 5.446753   formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q7 + Q8 + Q1:Q2 + Q3:Q4 + Q4:Q5 + Q5:Q6 + Q7:Q8 
##  ---Refining positions
##  ---  Moved a bit
##  ---Dropping Q2 
##     no.qtl =  7   pLOD = 5.014261   formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q7 + Q2:Q3 + Q3:Q4 + Q4:Q5 + Q6:Q7 
##  ---Refining positions
##  ---  Moved a bit
##  ---Dropping Q7 
##     no.qtl =  6   pLOD = 5.480166   formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q2:Q3 + Q3:Q4 + Q4:Q5 
##  ---Refining positions
##  ---  Moved a bit
##  ---Dropping Q5 
##     no.qtl =  5   pLOD = 6.575859   formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q2:Q3 + Q3:Q4 
##  ---Refining positions
##  ---  Moved a bit
## ** new best ** (pLOD increased by 0.1535)
##  ---Dropping Q5 
##     no.qtl =  4   pLOD = 6.830819   formula: y ~ Q1 + Q2 + Q3 + Q4 + Q2:Q3 + Q3:Q4 
##  ---Refining positions
##  ---  Moved a bit
## ** new best ** (pLOD increased by 0.0859)
##  ---Dropping Q2:Q3 
##     no.qtl =  4   pLOD = 6.603798   formula: y ~ Q1 + Q2 + Q3 + Q4 + Q3:Q4 
##  ---Refining positions
##  ---  Moved a bit
##  ---Dropping Q2 
##     no.qtl =  3   pLOD = 6.570747   formula: y ~ Q1 + Q2 + Q3 + Q2:Q3 
##  ---Refining positions
##  ---  Moved a bit
##  ---Dropping Q2 
##     no.qtl =  2   pLOD = 4.154861   formula: y ~ Q1 + Q2 
##  ---Refining positions
##  ---  Moved a bit
##  ---Dropping Q2 
##     no.qtl =  1   pLOD = 5.54997   formula: y ~ Q1 
##  ---Refining positions
##  ---  Moved a bit
##  ---One last pass through refineqtl

##   QTL object containing genotype probabilities. 
## 
##        name chr    pos n.gen
## Q1  2@107.0   2 107.00     3
## Q2  4@216.0   4 216.00     3
## Q3  12@68.0  12  68.00     3
## Q4 14@112.3  14 112.33     3
## 
##   Formula: y ~ Q1 + Q2 + Q3 + Q4 + Q2:Q3 + Q1:Q2 
## 
##   pLOD:  6.84

Four (4) covariates were defined.

Composite Interval Mapping (CIM) based on the number of covariates calculated above

out.cim <- cim(mimulus.qtl_cim, pheno.col=16, n.marcovar= 4, window=15, method="hk", error.prob=0.001, map.function="kosambi")
summary(out.cim)

##            chr   pos   lod
## c1.loc59     1  59.0 1.577
## c2.loc132    2 132.0 3.060
## c3.loc67     3  67.0 2.577
## c4.loc105    4 105.0 1.713
## AA268        5  90.1 1.624
## MgSTS545     6  88.4 1.050
## AAT296       7 164.1 1.446
## c8.loc91     8  91.0 1.555
## BD143        9 139.6 0.834
## c10.loc98   10  98.0 1.376
## c11.loc161  11 161.0 1.772
## c12.loc68   12  68.0 2.554
## c13.loc121  13 121.0 0.352
## AAT374      14 112.3 9.354

Permutation test for significant LOD threshold definition (using HK regression)

permutation_hk_cim<- cim(mimulus.qtl_cim, pheno.col=16, n.marcovar= 4, window=15, method="hk", error.prob=0.001, map.function="kosambi",n.perm=1000)
threshold_hk_cim<-summary(permutation_hk_cim)
summary(permutation_hk_cim)

## LOD thresholds (1000 permutations)
##     [,1]
## 5%  4.66
## 10% 4.21

Permutation distribution plot: significant threshold in percentile 95

plot.hk_cim<-plot(permutation_hk_cim, col="green")
abline(plot.hk_cim, v=threshold_hk_cim[1,1], col="red", lty=3, lwd=3)

Plot: Composite Interval Mapping (Haley-Knott method)
- Linkage Groups 1 to 7

plot(out.cim, chr=c(1,2,3,4,5,6,7), lty=1, col="green")
add.cim.covar(out.cim, chr=c(1,2,3,4,5,6,7)) 
abline(out.cim, h=threshold_hk_cim[1,1], col="red", lty=3)

Linkage Groups 8 to 14

plot(out.cim, chr=c(8,9,10,11,12,13,14), lty=1, col="green")
add.cim.covar(out.cim, chr=c(8,9,10,11,12,13,14)) 
abline(out.cim, h=threshold_hk_cim[1,1], col="red", lty=3)

QTL detection

summary_hk_cim <- summary(out.cim, perms=permutation_hk_cim, alpha = 0.05)
summary_hk_cim

##        chr pos  lod
## AAT374  14 112 9.35

QTL observation on the linkage map

qtl_cim <- makeqtl(mimulus.qtl_cim, chr=summary_hk_cim$chr, pos=summary_hk_cim$pos, what=c("draws","prob"))
plot(qtl_cim)

QTL additive and dominance effects detection

\(a=\displaystyle\frac{(\mu_{BB} - \mu_{AA})}{2}\)

\(d\)=\(\mu_{AB} - \displaystyle\frac{(\mu_{AA} + \mu_{BB})}{2}\)

Linkage Groups 1 to 7

chr1_7_cim <- effectscan(mimulus.qtl_cim, chr=c(1,2,3,4,5,6,7), pheno.col=16, get.se=FALSE, draw=TRUE,
gap=25, mtick=c("line","triangle"),add.legend=TRUE, alternate.chrid=FALSE, ylab="Effect", xlab="Linkage Group")

summary(chr1_7_cim)

##          chr pos     a       d
## c1.loc82   1  82 0.829  0.4676
## c2.loc43   2  43 0.985 -0.1917
## c3.loc89   3  89 1.233 -0.7740
## c4.loc36   4  36 0.510 -0.0782
## c5.loc90   5  90 0.632  0.4655
## c6.loc2    6   2 0.612 -0.6119
## c7.loc76   7  76 0.556  0.1220

+ **Linkage Groups 8 to 14**

chr8_14_cim <- effectscan(mimulus.qtl_cim, chr=c(8,9,10,11,12,13,14), pheno.col=16, get.se=FALSE, draw=TRUE,
gap=25, mtick=c("line","triangle"),add.legend=TRUE, alternate.chrid=FALSE, ylab="Effect", xlab="Linkage Group")

summary(chr8_14_cim)

##            chr pos     a       d
## c8.loc95     8  95 0.666 -0.1181
## c9.loc148    9 148 0.394  0.1512
## MgSTS622    10 133 0.672  0.4155
## MgSTS494    11 196 2.227 -1.6962
## c12.loc108  12 108 1.920  1.7240
## c13.loc120  13 120 0.387 -0.1104
## c14.loc114  14 114 1.865  0.0887

Position/marker presenting the highest effect in linkage group 11 (“MgSTS494”)

effectplot(f2_qtlmimulus, pheno.col=16, "MgSTS494")

Position/marker presenting the highest effect in linkage group 14 (“c14.loc114”)

effectplot(mimulus1, pheno.col=16, "c14.loc114")

Fitting the multiple-QTL model

out.fq_cim <- fitqtl(mimulus.qtl_cim, pheno.col=16, qtl=qtl_cim, method="imp", get.ests=TRUE)

## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 11 individuals with missing phenotypes.

summary(out.fq_cim)

## 
##      fitqtl summary
## 
## Method: multiple imputation 
## Model:  normal phenotype
## Number of observations : 276 
## 
## Full model result
## ----------------------------------  
## Model formula: y ~ Q1 
## 
##        df        SS         MS      LOD     %var Pvalue(Chi2)    Pvalue(F)
## Model   2  424.5818 212.290906 9.029189 13.98561 9.349995e-10 1.172073e-09
## Error 273 2611.2651   9.565074                                            
## Total 275 3035.8469                                                       
## 
## 
## Estimated effects:
## -----------------
##               est      SE      t
## Intercept 15.6535  0.1872 83.636
## 14@112.3a  1.8378  0.2759  6.662
## 14@112.3d  0.1112  0.3748  0.297

Table 2: Linkage Group, position, LOD score, QTLs effects (additive and dominant) and determination coefficient (%) for the Composite Interval Mapping (CIM) for the phenotype corolla width, using the Haley-Knott regression.

r2<-c(summary(out.fq_cim)$result.full[1,5])
rr<-as.numeric(r2)
table2<-data.frame("Linkage Group"=c(summary_hk_cim$chr[1]),
                   "Position"=c(summary_hk_cim$pos[1]),
                   "LOD"=c(round(summary_hk_cim$lod[1],4)),
                   "Additive effect"=c(round(summary(chr8_14_cim)$a[7],4)),
                   "Dominance Effect"=c(round(summary(chr8_14_cim)$d[7],3)),
                   "R2"=r2)
knitr::kable(table2)

Linkage.Group	Position	LOD	Additive.effect	Dominance.Effect	R2
14	112.3271	9.3542	1.865	0.089	13.98561

Conclusion (Composite Interval Mapping-CIM): Using the Haley-Knott method and considering the LOD threshold obtained through permutations, we observed 1 possible QTL for the phenotype corolla width in the linkage group 14, presenting mainly additive effects. As we can see by the determination coefficient, the percentage of explained phenotypic variance for the trait by the QTL was 13.64.

Comparison between Interval Mapping (IM) and Composite Interval Mapping (CIM)

Entire genome

plot(out_emmimulus1,out_hkmimulus1, out.cim,col=c("blue","red","green"))
abline(out.cim, h=threshold_hk_cim[1,1], col="dark green", lty=3)
abline(plot.hk.em, h= threshold_hk_im[1,1], col="red", lty=3)

Linkage Groups 1 to 7

plot(out_emmimulus1,out_hkmimulus1,chr=c(1,2,3,4,5,6,7), out.cim,col=c("blue", "red","green"))
abline(out.cim, h=threshold_hk_cim[1,1], col="dark green", lty=3)
abline(plot.hk.em, h= threshold_hk_im[1,1], col="red", lty=3)

Linkage Groups 8 to 14

plot(out_emmimulus1,out_hkmimulus1,chr=c(8,9,10,11,12,13,14), out.cim,col=c("blue","red","green"))
abline(out.cim, h=threshold_hk_cim[1,1], col="dark green", lty=3)
abline(plot.hk.em, h= threshold_hk_im[1,1], col="red", lty=3)

Linkage Group 14

plot(out_emmimulus1, out_hkmimulus1,chr=14,out.cim,col=c("blue","red","green"))
abline(out.cim, h=threshold_hk_cim[1,1], col="dark green", lty=3)
abline(plot.hk.em, h= threshold_hk_im[1,1], col= "red", lty=3)

Interval Mapping (IM) and Composite Interval Mapping (CIM) for QTL detection in a F2 Mimulus population

Alline Sekiya, Lillian Bibiano, Mariana Niederheitmann, Willian Giordani

June 7th, 2017

Interval Mapping (IM) for the phenotype corolla width

`plotMap(f2_qtlmimulus,main=substitute(paste("Genetic linkage map - F2 ")(italic('Mimulus guttatus x Mimulus nasutus'))))`

`plotMissing(f2_qtlmimulus, main="Missing genotypes")`

Corolla width (mm)

`#plot.pheno<-plot.pheno(f2_qtlmimulus, pheno.col = 16) plot(plot.pheno, col="yellow", xlab="Corola width (mm)", main="ww")`

Composite Interval Mapping (CIM) for the phenotype corolla width

Comparison between Interval Mapping (IM) and Composite Interval Mapping (CIM)