QTL analysis
Thomas Bergmann
23 September, 2021
1 Data preparation
Load libraries:
## Clear workspace
rm(list=ls())
## Libraries
library(qtl)
library(ggrepel)
library(dplyr)
library(ggplot2)
library(cowplot)
library(gdata)## gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.## ## gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.##
## Attaching package: 'gdata'## The following object is masked from 'package:gridExtra':
##
## combine## The following objects are masked from 'package:dplyr':
##
## combine, first, last## The following object is masked from 'package:stats':
##
## nobs## The following object is masked from 'package:utils':
##
## object.size## The following object is masked from 'package:base':
##
## startsWithlibrary(readxl)
library(gridExtra)
library(ggpubr)
library(VennDiagram)
library(grid)
library(futile.logger)
library(gplots)## Registered S3 method overwritten by 'gplots':
## method from
## reorder.factor gdata##
## Attaching package: 'gplots'## The following object is masked from 'package:gdata':
##
## reorder.factor## The following object is masked from 'package:stats':
##
## lowess## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyserPopulation structure:
- Two separate mapping populations: Population A / Population B
The genetic map was constructed according to Bromann via R/qtl. The QTL analysis is likewise done via R/qtl.
Data import:
## Load data
Population <- qtl::read.cross("csv",
"Data/",
"Population_Final.csv",
genotypes=c("AA","AB","BB"),
estimate.map=FALSE)## --Read the following data:
## 361 individuals
## 1118 markers
## 9 phenotypes
## --Cross type: f21.1 Quality check
Data check:
## F2 intercross
##
## No. individuals: 361
##
## No. phenotypes: 9
## Percent phenotyped: 100 46 46 50.1 43.5 47.4 50.1 44.3 46.8
##
## No. chromosomes: 10
## Autosomes: C01 C02 C03 C04a C04b C05 C06 C07 C08 C09
##
## Total markers: 1118
## No. markers: 121 122 233 29 36 131 131 142 139 34
## Percent genotyped: 99.8
## Genotypes (%): AA:26.4 AB:50.7 BB:23.0 not BB:0.0 not AA:0.0Genetic map:
## Make table
DT::datatable(GeneticMap <- as.data.frame(qtl::summaryMap(Population))) %>%
DT::formatRound(c("length","ave.spacing","max.spacing"), 3)Plot:
Note:
- LeafA/PetioleA –> Phenotypes in Population A
- LeafB1/B2/B3 and PetioleB1/B2/B3 –> Phenotypes in Population B
2 One-dimensional scan
Calculate genotype probabilities and load permutations:
## Calculate genotype probabilities
Population <- qtl::calc.genoprob(Population, step=1, error.prob=0.001)
## Simulate genotypes (almost no missing data here)
Population <- qtl::sim.geno(Population, n.draws=64, error.prob=0.001)
## Calculate significance threshold for phenotypes in population A
#operm.leafA.hk <- qtl::scanone(Population, pheno.col=2, method="hk", n.perm=10000)
#operm.petioleA.hk <- qtl::scanone(Population, pheno.col=3, method="hk", n.perm=10000)
## Calculate significance threshold for phenotypes in population B
#operm.leafB1.hk <- qtl::scanone(Population, pheno.col=4, method="hk", n.perm=10000)
#operm.leafB2.hk <- qtl::scanone(Population, pheno.col=5, method="hk", n.perm=10000)
#operm.leafB3.hk <- qtl::scanone(Population, pheno.col=6, method="hk", n.perm=10000)
#operm.petioleB1.hk <- qtl::scanone(Population, pheno.col=7, method="hk", n.perm=10000)
#operm.petioleB2.hk <- qtl::scanone(Population, pheno.col=8, method="hk", n.perm=10000)
#operm.petioleB3.hk <- qtl::scanone(Population, pheno.col=9, method="hk", n.perm=10000)
load("./Data/operm.leafA.hk.RData")
load("./Data/operm.petioleA.hk.RData")
load("./Data/operm.leafB1.hk.RData")
load("./Data/operm.leafB2.hk.RData")
load("./Data/operm.leafB3.hk.RData")
load("./Data/operm.petioleB1.hk.RData")
load("./Data/operm.petioleB2.hk.RData")
load("./Data/operm.petioleB3.hk.RData")Scanone:
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 195 individuals with missing phenotypes.## Petiole-scan in Population A
out.petioleA.hk <- qtl::scanone(Population, pheno.col=3, method="hk")## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 195 individuals with missing phenotypes.## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 180 individuals with missing phenotypes.## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 204 individuals with missing phenotypes.## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 190 individuals with missing phenotypes.## Petiole-scan1 in Population B
out.petioleB1.hk <- qtl::scanone(Population, pheno.col=7, method="hk")## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 180 individuals with missing phenotypes.## Petiole-scan2 in Population B
out.petioleB2.hk <- qtl::scanone(Population, pheno.col=8, method="hk")## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 201 individuals with missing phenotypes.## Petiole-scan3 in Population B
out.petioleB3.hk <- qtl::scanone(Population, pheno.col=9, method="hk")## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 192 individuals with missing phenotypes.2.1 Population A
2.1.1 Leaf-assay
## LOD thresholds (10000 permutations)
## lod
## 5% 3.61## There were no LOD peaks above the threshold.Non-parametric scan:
## Permutation
#operm.leafA.np <- qtl::scanone(Population, pheno.col=2, model="np", n.perm=10000)
load("./Data/operm.leafA.np.RData")
## Scan
out.leafA.np <- qtl::scanone(Population, pheno.col=2, model="np")## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 195 individuals with missing phenotypes.## LOD thresholds (10000 permutations)
## lod
## 5% 3.49## There were no LOD peaks above the threshold.LOD-profile:
## LOD-profiles
plot(out.leafA.hk, ylab="LOD score", ylim=c(0,7), col="deepskyblue1")
plot(out.leafA.np, ylab="LOD score", ylim=c(0,7), col="red", add=TRUE)
abline(h=3.61, col="deepskyblue1")
abline(h=3.49, col="firebrick")
LOD-profiles. Blue: Haley-Knott regression; red: non-parametric model
Note: Non-significant peak on C08.
2.1.2 Petiole-assay
## LOD thresholds (10000 permutations)
## lod
## 5% 3.65## chr pos lod pval
## Bn-scaff_19564_1-p17934 C01 53 6.23 3e-04LOD-profile:
## LOD-profile
plot(out.petioleA.hk, ylab="LOD score", ylim=c(0,7), col="deepskyblue1")
abline(h=3.65, col="deepskyblue1")
LOD-profile. Blue: Haley-Knott regression
QTL-scan controlling for the locus on linkage group C01:
## Get the peak marker on C01
mar <- qtl::find.marker(Population, "C01", 53)
g <- qtl::pull.geno(Population)[,mar]
## Check if there is missing data for the genotypes
sum(is.na(g))## [1] 0## Convert that it can be used as covariate
g <- cbind(as.numeric(g==1), as.numeric(g==2))
## Permutation
#operm.petioleA.hk.addcov <- scanone(Population, addcovar=g, pheno.col=3, chr="-C01", method="hk", n.perm=1000)
load("./Data/operm.petioleA.hk.addcov.RData")
## Re-scan
out.petioleA.hk.addcov <- scanone(Population, pheno.col=3, method="hk", addcovar=g)## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 195 individuals with missing phenotypes.## LOD thresholds (1000 permutations)
## lod
## 5% 3.61
## 10% 3.29## There were no LOD peaks above the threshold.QTL-scan controlling for the locus on linkage group C01 with interaction:
## Permutation
#operm.petioleA.hk.intcov <- scanone(Population, method="hk", pheno.col=3, n.perm=1000, addcovar=g, intcovar=g)
load("./Data/operm.petioleA.hk.intcov.RData")
## Re-scan
out.petioleA.hk.intcov <- scanone(Population, method="hk", pheno.col=3, addcovar=g, intcovar=g)## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 195 individuals with missing phenotypes.## LOD thresholds (1000 permutations)
## lod
## 5% 6.09
## 10% 5.69## There were no LOD peaks above the threshold.LOD-profiles:
## LOD-profile
plot(out.petioleA.hk, ylab="LOD score", ylim=c(0,7), col="deepskyblue1")
plot(out.petioleA.hk.addcov, col="darkslategray", ylab="LOD score", add=TRUE)
plot(out.petioleA.hk.intcov, col="darkorchid", ylab="LOD score", add=TRUE)
abline(h=3.65, col="deepskyblue1")
abline(h=3.61, col="darkslategray")
abline(h=6.09, col="darkorchid")
LOD-profiles. Blue = Haley-Knott regression; gray = Haley-Knott regression (additive covariate); purple = Haley-Knott regression (interactive covariate)
2.2 Population B
2.2.1 First leaf-assay
## LOD thresholds (10000 permutations)
## lod
## 5% 3.63## There were no LOD peaks above the threshold.Non-parametric scan:
## Permutation
#operm.leafB1.np <- qtl::scanone(Population, pheno.col=4, model="np", n.perm=10000)
load("./Data/operm.leafB1.np.RData")
## Scan
out.leafB1.np <- qtl::scanone(Population, pheno.col=4, model="np")## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 180 individuals with missing phenotypes.## LOD thresholds (10000 permutations)
## lod
## 5% 3.47## There were no LOD peaks above the threshold.LOD-profiles:
## LOD-profiles
plot(out.leafB1.hk, ylab="LOD score", ylim=c(0,7), col="deepskyblue1")
plot(out.leafB1.np, ylab="LOD score", ylim=c(0,7), col="red", add=TRUE)
abline(h=3.63, col="deepskyblue1")
abline(h=3.47, col="firebrick")
LOD-profiles. Blue: Haley-Knott regression; red: non-parametric model
2.2.2 Second leaf-assay
## LOD thresholds (10000 permutations)
## lod
## 5% 3.64## chr pos lod pval
## cC01.loc51 C01 51 6.16 3e-04LOD-profile:
## LOD-profile
plot(out.leafB2.hk, ylab="LOD score", ylim=c(0,7), col="deepskyblue1")
abline(h=3.64, col="deepskyblue1")
LOD-profile. Blue: Haley-Knott regression
QTL-scan controlling for the locus on linkage group C01:
## Get the peak marker on C01
mar <- qtl::find.marker(Population, "C01", 51)
g <- qtl::pull.geno(Population)[,mar]
## Check if there is missing data for the genotypes
sum(is.na(g))## [1] 0## Convert that it can be used as covariate
g <- cbind(as.numeric(g==1), as.numeric(g==2))
## Permutation
#operm.leafB2.hk.addcov <- scanone(Population, addcovar=g, pheno.col=5, chr="-C01", method="hk", n.perm=1000)
load("./Data/operm.leafB2.hk.addcov.RData")
## Re-scan
out.leafB2.hk.addcov <- scanone(Population, pheno.col=5, method="hk", addcovar=g)## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 204 individuals with missing phenotypes.## LOD thresholds (1000 permutations)
## lod
## 5% 3.66
## 10% 3.24## There were no LOD peaks above the threshold.QTL-scan controlling for the locus on linkage group C01 with interaction:
## Permutation
#operm.leafB2.hk.intcov <- scanone(Population, method="hk", pheno.col=5, n.perm=1000, addcovar=g, intcovar=g)
load("./Data/operm.leafB2.hk.intcov.RData")
## Re-scan
out.leafB2.hk.intcov <- scanone(Population, method="hk", pheno.col=5, addcovar=g, intcovar=g)## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 204 individuals with missing phenotypes.## LOD thresholds (1000 permutations)
## lod
## 5% 5.50
## 10% 5.13## chr pos lod pval
## Bn-scaff_15838_1-p1295927 C01 5.56 5.38 0.066LOD-profiles:
## LOD-profiles
plot(out.leafB2.hk, ylab="LOD score", ylim=c(0,7), col="deepskyblue1")
plot(out.leafB2.hk.addcov, col="darkslategray", ylab="LOD score", add=TRUE)
plot(out.leafB2.hk.intcov, col="darkorchid", ylab="LOD score", add=TRUE)
abline(h=3.64, col="deepskyblue1")
abline(h=3.66, col="darkslategray")
abline(h=5.5, col="darkorchid")
LOD-profiles. Blue = Haley-Knott regression; gray = Haley-Knott regression (additive covariate); purple = Haley-Knott regression (interactive covariate)
Non-significant peaks:
## chr pos lod
## Bn-scaff_16565_1-p899839 C02 46.6 5.02## chr pos lod
## cC03.loc26 C03 26 3.772.2.3 Third leaf-assay
## LOD thresholds (10000 permutations)
## lod
## 5% 3.59## chr pos lod pval
## Bn-scaff_16110_1-p976517 C07 79.6 4.07 0.0182LOD-profile:
## LOD-profile
plot(out.leafB3.hk, ylab="LOD score", ylim=c(0,7), col="deepskyblue1")
abline(h=3.59)
LOD-profile. Blue: Haley-Knott regression
QTL-scan controlling for the locus on linkage group C07:
## Get the peak marker on C07
mar <- qtl::find.marker(Population, "C07", 79.6)
g <- qtl::pull.geno(Population)[,mar]
## Chekc if there is missing data for the genotypes
sum(is.na(g))## [1] 0## Convert that it can be used as covariate
g <- cbind(as.numeric(g==1), as.numeric(g==2))
## Permutation
#operm.leafB3.hk.addcov <- scanone(Population, addcovar=g, pheno.col=6, chr="-C07", method="hk", n.perm=1000)
load("./Data/operm.leafB3.hk.addcov.RData")
## Re-scan
out.leafB3.hk.addcov <- scanone(Population, pheno.col=6, method="hk", addcovar=g)## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 190 individuals with missing phenotypes.## LOD thresholds (1000 permutations)
## lod
## 5% 3.68
## 10% 3.25## There were no LOD peaks above the threshold.QTL-scan controlling for the locus on linkage group C07 with interaction:
## Permutation
#operm.leafB3.hk.intcov <- scanone(Population, method="hk", pheno.col=6, n.perm=1000, addcovar=g, intcovar=g)
load("./Data/operm.leafB3.hk.intcov.RData")
## Re-scan
out.leafB3.hk.intcov <- scanone(Population, method="hk", pheno.col=6, addcovar=g, intcovar=g)## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 190 individuals with missing phenotypes.## LOD thresholds (1000 permutations)
## lod
## 5% 5.96
## 10% 5.50## There were no LOD peaks above the threshold.LOD-profiles:
## LOD-profiles
plot(out.leafB3.hk, ylab="LOD score", ylim=c(0,8), col="deepskyblue1")
plot(out.leafB3.hk.addcov, col="darkslategray", ylab="LOD score", add=TRUE)
plot(out.leafB3.hk.intcov, col="darkorchid", ylab="LOD score", add=TRUE)
abline(h=3.59, col="deepskyblue1")
abline(h=3.68, col="darkslategray")
abline(h=5.96, col="darkorchid")
LOD-profiles. Blue = Haley-Knott regression (Single-QTL); gray = Haley-Knott regression (additive covariate); purple = Haley-Knott regression (interactive covariate)
2.2.4 First petiole-assay
## LOD thresholds (10000 permutations)
## lod
## 5% 3.56## chr pos lod pval
## cC03.loc12 C03 12.0 3.79 0.0321
## Bn-scaff_16069_1-p4306874 C07 54.3 4.56 0.0060LOD-profile:
## LOD-profile
plot(out.petioleB1.hk, ylab="LOD score", ylim=c(0,7), col="deepskyblue1")
abline(h=3.56)
LOD-profile. Blue: Haley-Knott regression
QTL-scan controlling for the locus on linkage group C07:
## Get the peak marker on C07
mar <- qtl::find.marker(Population, "C07", 54.3)
g <- qtl::pull.geno(Population)[,mar]
## Check if there is missing datat for the genotypes
sum(is.na(g))## [1] 0## Convert that it can be used as covariate
g <- cbind(as.numeric(g==1), as.numeric(g==2))
## Permutation
#operm.petioleB1.hk.addcov <- scanone(Population, addcovar=g, pheno.col=7, chr="-C07", method="hk", n.perm=1000)
load("./Data/operm.petioleB1.hk.addcov.RData")
## Re-scan
out.petioleB1.hk.addcov <- scanone(Population, pheno.col=7, method="hk", addcovar=g)## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 180 individuals with missing phenotypes.## LOD thresholds (1000 permutations)
## lod
## 5% 3.61
## 10% 3.26## chr pos lod pval
## Bn-scaff_15877_1-p4822 C03 16.8 3.77 0.034
## Bn-scaff_16397_1-p165703 C06 47.8 3.50 0.062QTL-scan controlling for the locus on linkage group C07 with interaction:
## Re-scan
#operm.petioleB1.hk.intcov <- scanone(Population, method="hk", pheno.col=7, n.perm=1000, addcovar=g, intcovar=g)
load("./Data/operm.petioleB1.hk.intcov.RData")
out.petioleB1.hk.intcov <- scanone(Population, method="hk", pheno.col=7, addcovar=g, intcovar=g)## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 180 individuals with missing phenotypes.## LOD thresholds (1000 permutations)
## lod
## 5% 5.71
## 10% 5.35## chr pos lod pval
## Bn-scaff_15877_1-p4822 C03 16.8 5.76 0.042LOD-profiles:
## LOD-profiles
plot(out.petioleB1.hk, ylab="LOD score", ylim=c(0,8), col="deepskyblue1")
plot(out.petioleB1.hk.addcov, col="darkslategray", ylab="LOD score", add=TRUE)
plot(out.petioleB1.hk.intcov, col="darkorchid", ylab="LOD score", add=TRUE)
abline(h=3.56, col="deepskyblue1")
abline(h=3.61, col="darkslategray")
abline(h=5.71, col="darkorchid")
LOD-profiles. Blue = Haley-Knott regression; gray = Haley-Knott regression (additive covariate); purple = Haley-Knott regression (interactive covariate)
Non-significant peak:
## chr pos lod
## Bn-scaff_16397_1-p165703 C06 47.8 2.772.2.5 Second petiole-assay
## LOD thresholds (10000 permutations)
## lod
## 5% 3.66## There were no LOD peaks above the threshold.Non-parametric scan:
## Permutation
#operm.petioleB2.np <- qtl::scanone(Population, pheno.col=8, model="np", n.perm=1000)
load("./Data/operm.petioleB2.np.RData")
## Scan
out.petioleB2.np <- qtl::scanone(Population, pheno.col=8, model="np")## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 201 individuals with missing phenotypes.## LOD thresholds (1000 permutations)
## lod
## 5% 3.49## There were no LOD peaks above the threshold.LOD-profiles:
## LOD-profiles
plot(out.petioleB2.hk, ylab="LOD score", ylim=c(0,7), col="deepskyblue1")
plot(out.petioleB2.np, ylab="LOD score", ylim=c(0,7), col="red", add=TRUE)
abline(h=3.66, col="deepskyblue1")
abline(h=3.49, col="red")
LOD-profiles. Blue: Haley-Knott regression; red: non-parametric model
2.2.6 Third petiole-assay
## LOD thresholds (10000 permutations)
## lod
## 5% 3.58## There were no LOD peaks above the threshold.Non-parametric scan:
## Permutation
#operm.petioleB3.np <- qtl::scanone(Population, pheno.col=9, model="np", n.perm=1000)
load("./Data/operm.petioleB3.np.RData")
## Scan
out.petioleB3.np <- qtl::scanone(Population, pheno.col=9, model="np")## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 192 individuals with missing phenotypes.## LOD thresholds (1000 permutations)
## lod
## 5% 3.49## chr pos lod pval
## Bn-scaff_18936_1-p269153 C03 13.1 3.22 0.087
## Bn-scaff_16069_1-p2611780 C07 46.9 3.78 0.031LOD-profiles:
## LOD-profiles
plot(out.petioleB3.hk, ylab="LOD score", ylim=c(0,7), col="deepskyblue1")
plot(out.petioleB3.np, ylab="LOD score", ylim=c(0,7), col="firebrick", add=TRUE)
abline(h=3.58, col="deepskyblue1")
abline(h=3.49, col="firebrick")
LOD-profiles. Blue: Haley-Knott regression; red: non-parametric model
## Get the peak marker on C07
mar <- qtl::find.marker(Population, "C07", 46.9)
g <- qtl::pull.geno(Population)[,mar]
## Check if there is missing data for the genotypes
sum(is.na(g))## [1] 0## Convert that it can be used as covariate
g <- cbind(as.numeric(g==1), as.numeric(g==2))
## Re-scan
#operm.petioleB3.hk.addcov <- scanone(Population, addcovar=g, pheno.col=9, chr="-C07", method="hk", n.perm=1000)
load("./Data/operm.petioleB3.hk.addcov.RData")
out.petioleB3.hk.addcov <- scanone(Population, pheno.col=9, method="hk", addcovar=g)## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 192 individuals with missing phenotypes.## LOD thresholds (1000 permutations)
## lod
## 5% 3.52
## 10% 3.26## There were no LOD peaks above the threshold.## Re-scan
#operm.petioleB3.hk.intcov <- scanone(Population, method="hk", pheno.col=9, n.perm=1000, addcovar=g, intcovar=g)
load("./Data/operm.petioleB3.hk.intcov.RData")
out.petioleB3.hk.intcov <- scanone(Population, method="hk", pheno.col=9, addcovar=g, intcovar=g)## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 192 individuals with missing phenotypes.## LOD thresholds (1000 permutations)
## lod
## 5% 5.77
## 10% 5.41## There were no LOD peaks above the threshold.LOD-profiles:
## LOD-profiles
plot(out.petioleB3.hk, ylab="LOD score", ylim=c(0,8), col="deepskyblue1")
plot(out.petioleB3.np, col="firebrick", add=TRUE)
plot(out.petioleB3.hk.addcov, col="darkslategray", add=TRUE)
plot(out.petioleB3.hk.intcov, col="darkorchid", add=TRUE)
abline(h=3.58, col="deepskyblue1")
abline(h=3.49, col="firebrick")
abline(h=3.52, col="darkslategray")
abline(h=5.77, col="darkorchid")
LOD-profiles. Blue: Haley-Knott regression; red: non-parametric model; gray: Haley-Knott regression (additive model); purple: Haley-Knott regression (interactive model
2.3 Summary
LOD-profiles:
## Prepare window
par(mfrow=c(2,1))
## LOD-profiles leaf-assays
plot(out.leafA.hk, ylab="LOD score", col="red", ylim=c(0,7))
plot(out.leafB1.hk, ylab="LOD score", col="darkslategray", add=TRUE)
plot(out.leafB2.hk, ylab="LOD score", col="deepskyblue4", add=TRUE)
plot(out.leafB3.hk, ylab="LOD score", col="deepskyblue1", add=TRUE)
abline(h=3.5, col="gray30")
## LOD-profiles petiole-assays
plot(out.petioleA.hk, ylab="LOD score", ylim=c(0,7), col="red")
plot(out.petioleB1.hk, ylab="LOD score", col="darkslategray", add=TRUE)
plot(out.petioleB2.hk, ylab="LOD score", col="deepskyblue4", add=TRUE)
plot(out.petioleB3.np, ylab="LOD score", col="deepskyblue1", add=TRUE)
abline(h=3.5, col="gray30")
LOD-profiles. Top: leaf-assays. Bottom: petiole-assays. Red: LOD in Population A; blue: LOD in Population B. Blue shades indicate replications in Population B
3 Two-dimensional scan
3.1 Population A
3.1.1 Petiole-leaf assay
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 195 individuals with missing phenotypes.## --Running scanone
## --Running scantwo
## (C01,C01)
## (C01,C02)
## (C01,C03)
## (C01,C04a)
## (C01,C04b)
## (C01,C05)
## (C01,C06)
## (C01,C07)
## (C01,C08)
## (C01,C09)
## (C02,C02)
## (C02,C03)
## (C02,C04a)
## (C02,C04b)
## (C02,C05)
## (C02,C06)
## (C02,C07)
## (C02,C08)
## (C02,C09)
## (C03,C03)
## (C03,C04a)
## (C03,C04b)
## (C03,C05)
## (C03,C06)
## (C03,C07)
## (C03,C08)
## (C03,C09)
## (C04a,C04a)
## (C04a,C04b)
## (C04a,C05)
## (C04a,C06)
## (C04a,C07)
## (C04a,C08)
## (C04a,C09)
## (C04b,C04b)
## (C04b,C05)
## (C04b,C06)
## (C04b,C07)
## (C04b,C08)
## (C04b,C09)
## (C05,C05)
## (C05,C06)
## (C05,C07)
## (C05,C08)
## (C05,C09)
## (C06,C06)
## (C06,C07)
## (C06,C08)
## (C06,C09)
## (C07,C07)
## (C07,C08)
## (C07,C09)
## (C08,C08)
## (C08,C09)
## (C09,C09)load("./Data/operm2a.petioleA.hk.RData")
load("./Data/operm2b.petioleA.hk.RData")
operm2.petioleA.hk <- c(operm2a.petioleA.hk, operm2b.petioleA.hk)Result:
## Result
summary(out2.petioleA.hk, perms=operm2.petioleA.hk, alphas=c(0.2,0.2,0,0.2,0.2), pvalues=TRUE)## pos1f pos2f lod.full pval lod.fv1 pval lod.int pval pos1a pos2a
## cC02:cC03 51 49 7.77 0.294 4.26 0.998 1.17 1 67 49
## lod.add pval lod.av1 pval
## cC02:cC03 6.61 0.028 3.09 0.155No evidence for additional QTL.
3.2 Population B
3.2.1 Second-leaf assay
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 204 individuals with missing phenotypes.## --Running scanone
## --Running scantwo
## (C01,C01)
## (C01,C02)
## (C01,C03)
## (C01,C04a)
## (C01,C04b)
## (C01,C05)
## (C01,C06)
## (C01,C07)
## (C01,C08)
## (C01,C09)
## (C02,C02)
## (C02,C03)
## (C02,C04a)
## (C02,C04b)
## (C02,C05)
## (C02,C06)
## (C02,C07)
## (C02,C08)
## (C02,C09)
## (C03,C03)
## (C03,C04a)
## (C03,C04b)
## (C03,C05)
## (C03,C06)
## (C03,C07)
## (C03,C08)
## (C03,C09)
## (C04a,C04a)
## (C04a,C04b)
## (C04a,C05)
## (C04a,C06)
## (C04a,C07)
## (C04a,C08)
## (C04a,C09)
## (C04b,C04b)
## (C04b,C05)
## (C04b,C06)
## (C04b,C07)
## (C04b,C08)
## (C04b,C09)
## (C05,C05)
## (C05,C06)
## (C05,C07)
## (C05,C08)
## (C05,C09)
## (C06,C06)
## (C06,C07)
## (C06,C08)
## (C06,C09)
## (C07,C07)
## (C07,C08)
## (C07,C09)
## (C08,C08)
## (C08,C09)
## (C09,C09)## Permutation
load("./Data/operm2a.leafB2.hk.RData")
load("./Data/operm2b.leafB2.hk.RData")
operm2.leafB2.hk <- c(operm2a.leafB2.hk, operm2b.leafB2.hk)Result:
## Result
summary(out2.leafB2.hk, perms=operm2.leafB2.hk, alphas=c(0.2,0.2,0,0.2,0.2), pvalues=TRUE)## There were no pairs of loci meeting the criteria.3.2.2 Third-leaf assay
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 190 individuals with missing phenotypes.## --Running scanone
## --Running scantwo
## (C01,C01)
## (C01,C02)
## (C01,C03)
## (C01,C04a)
## (C01,C04b)
## (C01,C05)
## (C01,C06)
## (C01,C07)
## (C01,C08)
## (C01,C09)
## (C02,C02)
## (C02,C03)
## (C02,C04a)
## (C02,C04b)
## (C02,C05)
## (C02,C06)
## (C02,C07)
## (C02,C08)
## (C02,C09)
## (C03,C03)
## (C03,C04a)
## (C03,C04b)
## (C03,C05)
## (C03,C06)
## (C03,C07)
## (C03,C08)
## (C03,C09)
## (C04a,C04a)
## (C04a,C04b)
## (C04a,C05)
## (C04a,C06)
## (C04a,C07)
## (C04a,C08)
## (C04a,C09)
## (C04b,C04b)
## (C04b,C05)
## (C04b,C06)
## (C04b,C07)
## (C04b,C08)
## (C04b,C09)
## (C05,C05)
## (C05,C06)
## (C05,C07)
## (C05,C08)
## (C05,C09)
## (C06,C06)
## (C06,C07)
## (C06,C08)
## (C06,C09)
## (C07,C07)
## (C07,C08)
## (C07,C09)
## (C08,C08)
## (C08,C09)
## (C09,C09)## Permutation
load("./Data/operm2a.leafB3.hk.RData")
load("./Data/operm2b.leafB3.hk.RData")
operm2.leafB3.hk <- c(operm2a.leafB3.hk, operm2b.leafB3.hk)Result:
## Result
summary(out2.leafB3.hk, perms=operm2.leafB3.hk, alphas=c(0.2,0.2,0,0.2,0.2), pvalues=TRUE)## There were no pairs of loci meeting the criteria.3.2.3 First-petiole assay
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 180 individuals with missing phenotypes.## --Running scanone
## --Running scantwo
## (C01,C01)
## (C01,C02)
## (C01,C03)
## (C01,C04a)
## (C01,C04b)
## (C01,C05)
## (C01,C06)
## (C01,C07)
## (C01,C08)
## (C01,C09)
## (C02,C02)
## (C02,C03)
## (C02,C04a)
## (C02,C04b)
## (C02,C05)
## (C02,C06)
## (C02,C07)
## (C02,C08)
## (C02,C09)
## (C03,C03)
## (C03,C04a)
## (C03,C04b)
## (C03,C05)
## (C03,C06)
## (C03,C07)
## (C03,C08)
## (C03,C09)
## (C04a,C04a)
## (C04a,C04b)
## (C04a,C05)
## (C04a,C06)
## (C04a,C07)
## (C04a,C08)
## (C04a,C09)
## (C04b,C04b)
## (C04b,C05)
## (C04b,C06)
## (C04b,C07)
## (C04b,C08)
## (C04b,C09)
## (C05,C05)
## (C05,C06)
## (C05,C07)
## (C05,C08)
## (C05,C09)
## (C06,C06)
## (C06,C07)
## (C06,C08)
## (C06,C09)
## (C07,C07)
## (C07,C08)
## (C07,C09)
## (C08,C08)
## (C08,C09)
## (C09,C09)## Permutation
load("./Data/operm2a.petioleB1.hk.RData")
load("./Data/operm2b.petioleB1.hk.RData")
operm2.petioleB1.hk <- c(operm2a.petioleB1.hk, operm2b.petioleB1.hk)Result:
## Result
summary(out2.petioleB1.hk, perms=operm2.petioleB1.hk, alphas=c(0.2,0.2,0,0.2,0.2), pvalues=TRUE)## pos1f pos2f lod.full pval lod.fv1 pval lod.int pval pos1a pos2a
## cC03:cC07 0 54 9.61 0.012 5.43 0.578 2.036 1 17 54
## cC06:cC07 45 54 8.54 0.102 4.36 0.998 0.914 1 48 54
## lod.add pval lod.av1 pval
## cC03:cC07 7.57 0.0075 3.39 0.124
## cC06:cC07 7.62 0.0065 3.44 0.114Evidence for the QTL on C06 is weak.
3.2.4 Third-petiole assay
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 192 individuals with missing phenotypes.## --Running scanone
## --Running scantwo
## (C01,C01)
## (C01,C02)
## (C01,C03)
## (C01,C04a)
## (C01,C04b)
## (C01,C05)
## (C01,C06)
## (C01,C07)
## (C01,C08)
## (C01,C09)
## (C02,C02)
## (C02,C03)
## (C02,C04a)
## (C02,C04b)
## (C02,C05)
## (C02,C06)
## (C02,C07)
## (C02,C08)
## (C02,C09)
## (C03,C03)
## (C03,C04a)
## (C03,C04b)
## (C03,C05)
## (C03,C06)
## (C03,C07)
## (C03,C08)
## (C03,C09)
## (C04a,C04a)
## (C04a,C04b)
## (C04a,C05)
## (C04a,C06)
## (C04a,C07)
## (C04a,C08)
## (C04a,C09)
## (C04b,C04b)
## (C04b,C05)
## (C04b,C06)
## (C04b,C07)
## (C04b,C08)
## (C04b,C09)
## (C05,C05)
## (C05,C06)
## (C05,C07)
## (C05,C08)
## (C05,C09)
## (C06,C06)
## (C06,C07)
## (C06,C08)
## (C06,C09)
## (C07,C07)
## (C07,C08)
## (C07,C09)
## (C08,C08)
## (C08,C09)
## (C09,C09)## Permutation
load("./Data/operm2a.petioleB3.hk.RData")
load("./Data/operm2b.petioleB3.hk.RData")
operm2.petioleB3.hk <- c(operm2a.petioleB3.hk, operm2b.petioleB3.hk)Result:
## Result
summary(out2.petioleB3.hk, perms=operm2b.petioleB3.hk, alphas=c(0.4,0.4,0,0.4,0.4), pvalues=TRUE)## pos1f pos2f lod.full pval lod.fv1 pval lod.int pval pos1a pos2a
## cC03:cC07 70 63 7.55 0.444 4.70 0.974 1.84 1 91 59
## cC06:cC07 47 47 7.48 0.482 4.63 0.979 1.52 1 47 58
## cC07:cC07 38 45 6.82 0.794 3.98 1.000 0.24 1 40 44
## lod.add pval lod.av1 pval
## cC03:cC07 5.71 0.119 2.86 0.365
## cC06:cC07 5.96 0.086 3.12 0.219
## cC07:cC07 6.58 0.033 3.74 0.053Though not significant due to the different method (Haley-knott instead of the non-parametric method), the locus on C06 appears again.
3.3 Summary
A weak interaction between the QTL on C03 and C07 was identified:
3.3.1 C03:C07-Interaction
First petiole-assay
## PetioleB1
PhenosPetioleB1 <- qtl::pull.pheno(Population, pheno.col=7)
PhenosPetioleB1 <- PhenosPetioleB1[181:361]
## Get markers
mar <- qtl::find.marker(Population, chr=c("C03", "C07"), c(12,54))
par(mfrow=c(1,2))
qtl::plotPXG(Population, marker=mar, 7)
qtl::effectplot(Population, pheno.col=7, main="Interaction C03:C07", mname1="C03@12", mname2="C07@54", ylim=c(21,35))
Effectplot for the QTL on C03 and C07; C03 = Bn-scaff_18936_1-p269153, C07 = Bn-scaff_16069_1-p4252709
The plot shows the weak effect of the QTL on C03. There are not enough individuals for AA:BB and BB:AA genotypes. The effect is not clear. We assume rather additivity of these two QTL.
Note:
BRA1896: A-allele
BRA1909: B-allele
Third-petiole-assay
## Get markers
mar <- qtl::find.marker(Population, chr=c("C03", "C07"), c(13.1,46.9))
par(mfrow=c(1,2))
qtl::plotPXG(Population, marker=mar, 7)
qtl::effectplot(Population, pheno.col=7, main="Interaction C03:C07", mname1="C03@13.1", mname2="C07@46.9", ylim=c(24,35))
Effectplot for the QTL on C03 and C07; C03 = Bn-scaff_18936_1-p269153, C07 = Bn-scaff_16069_1-p2611780
Same pattern as before. The AA:BB and BB:AA genotypes are very rare. The effect is not clearly to determine. We assume additivity.
3.3.2 C06:C07-Interaction
## Get markers
mar <- qtl::find.marker(Population, chr=c("C07","C06"), c(54,47))
par(mfrow=c(1,2))
qtl::plotPXG(Population, marker=mar, 7)
qtl::effectplot(Population, pheno.col=7, main="Interaction C06:C07", mname1="C06@46.9", mname2="C07@54", ylim=c(20,35))
Effectplot for the QTL on C06 and C07; C06 = Bn-scaff_16397_1-p360292, C07 = Bn-scaff_16069_1-p4252709
We can see an additive effect of C06.
4 Multiple-QTL analysis
4.1 Population A
4.1.1 Petiole-assay
Define the initial QTL:
## QTL object containing genotype probabilities.
##
## name chr pos n.gen
## Q1 C01@53.0 C01 53 3Fit the model:
## fit the model
petioleA.fq <- qtl::fitqtl(Population, pheno.col=3, qtl=C01a, get.ests=TRUE, formula=y~Q1, method="hk")## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 195 individuals with missing phenotypes.##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 166
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1
##
## df SS MS LOD %var
## Model 2 1275.533 637.76671 6.215658 15.83868
## Error 163 6777.750 41.58129
## Total 165 8053.283
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 38.4778 0.5009 76.821
## C01@53.0a 3.8216 0.6959 5.491
## C01@53.0d -0.6285 1.0023 -0.627Check for additional/interacting QTL with the initial model:
## Scan for additional QTL and check
add.qtl <- qtl::addqtl(Population, qtl=C01a, method="hk", formula=y~Q1, pheno.col=3)## Warning in qtl::addqtl(Population, qtl = C01a, method = "hk", formula = y ~ : Dropping 195 individuals with missing phenotypes.## There were no LOD peaks above the threshold.There is no evidence for any additional QTL. We perform the last forward/backward stepwise approach to check for the best model.
Calculate penalties:
## main heavy light
## 3.567357 5.895156 3.314516Perform the search:
## Perform best model search
print(stepPetioleA <- qtl::stepwiseqtl(Population, qtl=C01a, max.qtl=6,
penalties=pen, verbose=TRUE, pheno.col=3,
method="hk", keeptrace=TRUE))## -Initial scan
## ---Starting at a model with 1 QTL
## ---Refining positions
## --- Moved a bit
## ** new best ** (pLOD increased by 2.6622)
## no.qtl = 1 pLOD = 2.662246 formula: y ~ Q1
## -Step 1
## ---Scanning for additive qtl
## plod = 2.001449
## ---Scanning for QTL interacting with Q1
## plod = 0.6854424
## ---Refining positions
## no.qtl = 2 pLOD = 2.001449 formula: y ~ Q1 + Q2
## -Step 2
## ---Scanning for additive qtl
## plod = 1.833052
## ---Scanning for QTL interacting with Q1
## plod = 0.4885735
## ---Scanning for QTL interacting with Q2
## plod = -0.698032
## ---Look for additional interactions
## plod = -0.1225414
## ---Refining positions
## no.qtl = 3 pLOD = 1.833052 formula: y ~ Q1 + Q2 + Q3
## -Step 3
## ---Scanning for additive qtl
## plod = 1.195204
## ---Scanning for QTL interacting with Q1
## plod = -0.8918062
## ---Scanning for QTL interacting with Q2
## plod = -1.139155
## ---Scanning for QTL interacting with Q3
## plod = -0.8426167
## ---Look for additional interactions
## plod = 0.02014208
## ---Refining positions
## --- Moved a bit
## no.qtl = 4 pLOD = 1.35658 formula: y ~ Q1 + Q2 + Q3 + Q4
## -Step 4
## ---Scanning for additive qtl
## plod = -0.01009704
## ---Scanning for QTL interacting with Q1
## plod = -0.2450709
## ---Scanning for QTL interacting with Q2
## plod = -2.183588
## ---Scanning for QTL interacting with Q3
## plod = -2.160282
## ---Scanning for QTL interacting with Q4
## plod = -1.359436
## ---Look for additional interactions
## plod = -0.7399734
## ---Refining positions
## --- Moved a bit
## no.qtl = 5 pLOD = 0.124667 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5
## -Step 5
## ---Scanning for additive qtl
## plod = -1.840412
## ---Scanning for QTL interacting with Q1
## plod = -0.8684455
## ---Scanning for QTL interacting with Q2
## plod = -3.148726
## ---Scanning for QTL interacting with Q3
## plod = -3.382375
## ---Scanning for QTL interacting with Q4
## plod = -3.271828
## ---Scanning for QTL interacting with Q5
## plod = -3.43732
## ---Look for additional interactions
## plod = -1.94628
## ---Refining positions
## --- Moved a bit
## no.qtl = 6 pLOD = -0.725889 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q1:Q6
## -Starting backward deletion
## ---Dropping Q6
## no.qtl = 5 pLOD = -0.01497077 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5
## ---Refining positions
## --- Moved a bit
## ---Dropping Q5
## no.qtl = 4 pLOD = 1.249795 formula: y ~ Q1 + Q2 + Q3 + Q4
## ---Refining positions
## --- Moved a bit
## ---Dropping Q4
## no.qtl = 3 pLOD = 1.740818 formula: y ~ Q1 + Q2 + Q3
## ---Refining positions
## --- Moved a bit
## ---Dropping Q3
## no.qtl = 2 pLOD = 2.001449 formula: y ~ Q1 + Q2
## ---Refining positions
## ---Dropping Q2
## no.qtl = 1 pLOD = 2.662246 formula: y ~ Q1
## ---Refining positions
## ---One last pass through refineqtl
## QTL object containing genotype probabilities.
##
## name chr pos n.gen
## Q1 C01@53.0 C01 52.98 3
##
## Formula: y ~ Q1
##
## pLOD: 2.662Best fit is for the QTL on C01. We adjust the model and check the effect:
## Refine the position
refC01a <- qtl::refineqtl(Population, pheno.col=3, method="hk", qtl=C01a, keeplodprofile=TRUE)## pos: 53
## Iteration 1
## Q1 pos: 53 -> 52.97992
## LOD increase: 0.014
## all pos: 53 -> 52.97992
## LOD increase at this iteration: 0.014
## Iteration 2
## Q1 pos: 52.97992 -> 52.97992
## LOD increase: 0
## all pos: 52.97992 -> 52.97992
## LOD increase at this iteration: 0
## overall pos: 53 -> 52.97992
## LOD increase overall: 0.014Fit the model:
## Fit the model
petioleA.final <- qtl::fitqtl(Population, pheno.col=3, qtl=refC01a, get.ests=TRUE, formula=y~Q1, method="hk")## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 195 individuals with missing phenotypes.##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 166
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1
##
## df SS MS LOD %var
## Model 2 1278.155 639.0775 6.229603 15.87123
## Error 163 6775.128 41.5652
## Total 165 8053.283
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 38.4788 0.5008 76.837
## C01@53.0a 3.8247 0.6954 5.500
## C01@53.0d -0.6203 1.0016 -0.6194.2 Population B
4.2.1 Second-leaf-assay
Define initial model:
## QTL object containing genotype probabilities.
##
## name chr pos n.gen
## Q1 C01@51.0 C01 51 3Fit the model:
## Fit the model
leafB2.fq <- qtl::fitqtl(Population, pheno.col=5, qtl=C01b, get.ests=TRUE, formula=y~Q1, method="hk")## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 204 individuals with missing phenotypes.##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 157
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1
##
## df SS MS LOD %var
## Model 2 1068560 534279.98 6.160004 16.53035
## Error 154 5395673 35036.83
## Total 156 6464233
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 630.12 15.52 40.613
## C01@51.0a 123.92 24.78 5.002
## C01@51.0d 88.16 31.36 2.811Check for additional QTL and interactions with the initial model:
## Scan for additional QTL and check
add.qtl <- qtl::addqtl(Population, qtl=C01b, method="hk", formula=y~Q1, pheno.col=5)## Warning in qtl::addqtl(Population, qtl = C01b, method = "hk", formula = y ~ : Dropping 204 individuals with missing phenotypes.## There were no LOD peaks above the threshold.Calculate the penalties:
## main heavy light
## 3.621574 6.032504 3.383352Perform the search:
## Search best model
print(stepleafB2 <- qtl::stepwiseqtl(Population, qtl=C01b, max.qtl=6,
penalties=pen, verbose=TRUE, pheno.col=5,
method="hk", keeptrace=TRUE))## -Initial scan
## ---Starting at a model with 1 QTL
## ---Refining positions
## ** new best ** (pLOD increased by 2.5384)
## no.qtl = 1 pLOD = 2.53843 formula: y ~ Q1
## -Step 1
## ---Scanning for additive qtl
## plod = 1.477771
## ---Scanning for QTL interacting with Q1
## plod = 1.123363
## ---Refining positions
## no.qtl = 2 pLOD = 1.477771 formula: y ~ Q1 + Q2
## -Step 2
## ---Scanning for additive qtl
## plod = 0.6545759
## ---Scanning for QTL interacting with Q1
## plod = -0.7507815
## ---Scanning for QTL interacting with Q2
## plod = -1.536238
## ---Look for additional interactions
## plod = 0.258366
## ---Refining positions
## no.qtl = 3 pLOD = 0.6545759 formula: y ~ Q1 + Q2 + Q3
## -Step 3
## ---Scanning for additive qtl
## plod = -1.181135
## ---Scanning for QTL interacting with Q1
## plod = -1.772117
## ---Scanning for QTL interacting with Q2
## plod = -2.460891
## ---Scanning for QTL interacting with Q3
## plod = -3.258482
## ---Look for additional interactions
## plod = -0.7810706
## ---Refining positions
## --- Moved a bit
## no.qtl = 3 pLOD = -0.4898662 formula: y ~ Q1 + Q2 + Q3 + Q1:Q2
## -Step 4
## ---Scanning for additive qtl
## plod = -2.131738
## ---Scanning for QTL interacting with Q1
## plod = -5.756506
## ---Scanning for QTL interacting with Q2
## plod = -6.090944
## ---Scanning for QTL interacting with Q3
## plod = -4.513639
## ---Look for additional interactions
## plod = -5.111015
## ---Refining positions
## --- Moved a bit
## no.qtl = 4 pLOD = -1.965212 formula: y ~ Q1 + Q2 + Q3 + Q1:Q2 + Q4
## -Step 5
## ---Scanning for additive qtl
## plod = -3.632204
## ---Scanning for QTL interacting with Q1
## plod = -7.015901
## ---Scanning for QTL interacting with Q2
## plod = -7.494048
## ---Scanning for QTL interacting with Q3
## plod = -5.363201
## ---Scanning for QTL interacting with Q4
## plod = -4.630628
## ---Look for additional interactions
## plod = -6.990154
## ---Refining positions
## --- Moved a bit
## no.qtl = 5 pLOD = -3.568155 formula: y ~ Q1 + Q2 + Q3 + Q1:Q2 + Q4 + Q5
## -Step 6
## ---Scanning for additive qtl
## plod = -5.467949
## ---Scanning for QTL interacting with Q1
## plod = -8.955501
## ---Scanning for QTL interacting with Q2
## plod = -8.33717
## ---Scanning for QTL interacting with Q3
## plod = -7.203786
## ---Scanning for QTL interacting with Q4
## plod = -6.363483
## ---Scanning for QTL interacting with Q5
## plod = -6.411639
## ---Look for additional interactions
## plod = -5.967087
## ---Refining positions
## --- Moved a bit
## no.qtl = 6 pLOD = -5.291569 formula: y ~ Q1 + Q2 + Q3 + Q1:Q2 + Q4 + Q5 + Q6
## -Starting backward deletion
## ---Dropping Q2
## no.qtl = 5 pLOD = -3.064144 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5
## ---Refining positions
## --- Moved a bit
## ---Dropping Q3
## no.qtl = 4 pLOD = -1.448396 formula: y ~ Q1 + Q2 + Q3 + Q4
## ---Refining positions
## --- Moved a bit
## ---Dropping Q3
## no.qtl = 3 pLOD = -0.02659778 formula: y ~ Q1 + Q2 + Q3
## ---Refining positions
## --- Moved a bit
## ---Dropping Q2
## no.qtl = 2 pLOD = 1.477771 formula: y ~ Q1 + Q2
## ---Refining positions
## ---Dropping Q2
## no.qtl = 1 pLOD = 2.53843 formula: y ~ Q1
## ---Refining positions
## ---One last pass through refineqtl
## QTL object containing genotype probabilities.
##
## name chr pos n.gen
## Q1 C01@51.0 C01 51 3
##
## Formula: y ~ Q1
##
## pLOD: 2.538Best fit is for the QTL on C01. We adjust the model and check the effect:
## pos: 51
## Iteration 1
## Q1 pos: 51 -> 51
## LOD increase: 0
## all pos: 51 -> 51
## LOD increase at this iteration: 0
## overall pos: 51 -> 51
## LOD increase overall: 0Fit the model:
## Fit the model
leafB2.final <- qtl::fitqtl(Population, pheno.col=5, qtl=refC01b, get.ests=TRUE, formula=y~Q1, method="hk")## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 204 individuals with missing phenotypes.##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 157
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1
##
## df SS MS LOD %var
## Model 2 1068560 534279.98 6.160004 16.53035
## Error 154 5395673 35036.83
## Total 156 6464233
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 630.12 15.52 40.613
## C01@51.0a 123.92 24.78 5.002
## C01@51.0d 88.16 31.36 2.8114.2.2 Third-leaf-assay
Define initial model:
## Make the model
print(leafC07 <- qtl::makeqtl(Population, chr=c("C07"), pos=c(79.6), what="prob"))## QTL object containing genotype probabilities.
##
## name chr pos n.gen
## Q1 C07@79.6 C07 79.557 3Fit the model:
## Model-fit
leafB3.fq <- qtl::fitqtl(Population, pheno.col=6, qtl=leafC07, get.ests=TRUE, formula=y~Q1, method="hk")## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 190 individuals with missing phenotypes.##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 171
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1
##
## df SS MS LOD %var
## Model 2 378360.9 189180.4 4.072145 10.38668
## Error 168 3264390.8 19430.9
## Total 170 3642751.7
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 418.60 10.68 39.212
## C07@79.6a 15.92 15.25 1.044
## C07@79.6d -90.28 21.35 -4.229Check for additional QTL:
## Scan for additional QTL and check
add.qtl <- qtl::addqtl(Population, qtl=leafC07, method="hk", formula=y~Q1, pheno.col=6)## Warning in qtl::addqtl(Population, qtl = leafC07, method = "hk", formula = y ~ : Dropping 190 individuals with missing phenotypes.## There were no LOD peaks above the threshold.Calculate the penalties:
## main heavy light
## 3.556775 6.411282 4.542966Perform the search:
## Perform best model search
print(stepleafB3 <- qtl::stepwiseqtl(Population, qtl=leafC07, max.qtl=6,
penalties=pen, verbose=TRUE, pheno.col=6,
method="hk", keeptrace=TRUE))## -Initial scan
## ---Starting at a model with 1 QTL
## ---Refining positions
## --- Moved a bit
## ** new best ** (pLOD increased by 0.5154)
## no.qtl = 1 pLOD = 0.5153698 formula: y ~ Q1
## -Step 1
## ---Scanning for additive qtl
## plod = -1.210329
## ---Scanning for QTL interacting with Q1
## plod = -3.72423
## ---Refining positions
## no.qtl = 2 pLOD = -1.210329 formula: y ~ Q1 + Q2
## -Step 2
## ---Scanning for additive qtl
## plod = -3.12219
## ---Scanning for QTL interacting with Q1
## plod = -6.217956
## ---Scanning for QTL interacting with Q2
## plod = -3.296449
## ---Look for additional interactions
## plod = -3.72423
## ---Refining positions
## no.qtl = 3 pLOD = -3.12219 formula: y ~ Q1 + Q2 + Q3
## -Step 3
## ---Scanning for additive qtl
## plod = -5.083301
## ---Scanning for QTL interacting with Q1
## plod = -8.120255
## ---Scanning for QTL interacting with Q2
## plod = -5.214576
## ---Scanning for QTL interacting with Q3
## plod = -6.302031
## ---Look for additional interactions
## plod = -5.684215
## ---Refining positions
## no.qtl = 4 pLOD = -5.083301 formula: y ~ Q1 + Q2 + Q3 + Q4
## -Step 4
## ---Scanning for additive qtl
## plod = -6.929262
## ---Scanning for QTL interacting with Q1
## plod = -9.973298
## ---Scanning for QTL interacting with Q2
## plod = -7.312154
## ---Scanning for QTL interacting with Q3
## plod = -9.059819
## ---Scanning for QTL interacting with Q4
## plod = -8.826742
## ---Look for additional interactions
## plod = -6.302031
## ---Refining positions
## --- Moved a bit
## no.qtl = 4 pLOD = -6.270718 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q3:Q4
## -Step 5
## ---Scanning for additive qtl
## plod = -7.932176
## ---Scanning for QTL interacting with Q1
## plod = -11.14216
## ---Scanning for QTL interacting with Q2
## plod = -8.5479
## ---Scanning for QTL interacting with Q3
## plod = -12.79259
## ---Scanning for QTL interacting with Q4
## plod = -12.51797
## ---Look for additional interactions
## plod = -9.050759
## ---Refining positions
## --- Moved a bit
## no.qtl = 5 pLOD = -7.882919 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q3:Q4 + Q5
## -Step 6
## ---Scanning for additive qtl
## plod = -8.743833
## ---Scanning for QTL interacting with Q1
## plod = -12.80632
## ---Scanning for QTL interacting with Q2
## plod = -10.56033
## ---Scanning for QTL interacting with Q3
## plod = -14.06411
## ---Scanning for QTL interacting with Q4
## plod = -13.82161
## ---Scanning for QTL interacting with Q5
## plod = -10.80944
## ---Look for additional interactions
## plod = -12.27849
## ---Refining positions
## no.qtl = 6 pLOD = -8.743833 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q3:Q4 + Q5 + Q6
## -Starting backward deletion
## ---Dropping Q4
## no.qtl = 5 pLOD = -6.561247 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5
## ---Refining positions
## --- Moved a bit
## ---Dropping Q2
## no.qtl = 4 pLOD = -4.57595 formula: y ~ Q1 + Q2 + Q3 + Q4
## ---Refining positions
## --- Moved a bit
## ---Dropping Q2
## no.qtl = 3 pLOD = -3.040004 formula: y ~ Q1 + Q2 + Q3
## ---Refining positions
## --- Moved a bit
## ---Dropping Q3
## no.qtl = 2 pLOD = -2.945448 formula: y ~ Q1 + Q2
## ---Refining positions
## --- Moved a bit
## ---Dropping Q2
## no.qtl = 1 pLOD = 0.5153698 formula: y ~ Q1
## ---Refining positions
## --- Moved a bit
## ---One last pass through refineqtl
## QTL object containing genotype probabilities.
##
## name chr pos n.gen
## Q1 C07@79.6 C07 79.557 3
##
## Formula: y ~ Q1
##
## pLOD: 0.515We adjust the model and check the effect:
## refine the position
refleafC07 <- qtl::refineqtl(Population, pheno.col=6, method="hk", qtl=leafC07)## pos: 79.55689
## Iteration 1
## Q1 pos: 79.55689 -> 79.55689
## LOD increase: 0
## all pos: 79.55689 -> 79.55689
## LOD increase at this iteration: 0
## Iteration 2
## Q1 pos: 79.55689 -> 79.55689
## LOD increase: 0
## all pos: 79.55689 -> 79.55689
## LOD increase at this iteration: 0
## overall pos: 79.55689 -> 79.55689
## LOD increase overall: 0Fit the model:
## Model-fit
leafB3.final <- qtl::fitqtl(Population, pheno.col=6, qtl=refleafC07, get.ests=TRUE, formula=y~Q1, method="hk")## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 190 individuals with missing phenotypes.##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 171
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1
##
## df SS MS LOD %var
## Model 2 378360.9 189180.4 4.072145 10.38668
## Error 168 3264390.8 19430.9
## Total 170 3642751.7
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 418.60 10.68 39.212
## C07@79.6a 15.92 15.25 1.044
## C07@79.6d -90.28 21.35 -4.229The QTL explains about 10 % of variance. The additive ffect is 15.92 and the dominance effect comes from BRA1896 (resistant) and is -90.28. We take a look at the phenotype data:
## Get markers
mar <- qtl::find.marker(Population, chr=c("C07"), c(79.6))
par(mfrow=c(1,2))
qtl::plotPXG(Population, marker=mar, 6)
qtl::effectplot(Population, main="Effectplot", pheno.col=6, mname1="C07@79.6")
Effectplot for l3QTLb
Heterozygous genotypes have a higher level of resistance than the homozygous AA-individuals. AA-individuals are less than AB- and BB- genotypes. Measurement in AA may be inaccurate due to this smaller group and the variation in the leaf-assay.
4.2.3 First petiole-assay
Define initial model:
## Make the model
print(twoQTL <- qtl::makeqtl(Population, chr=c("C03","C07"), pos=c(12,54.2), what="prob"))## QTL object containing genotype probabilities.
##
## name chr pos n.gen
## Q1 C03@12.0 C03 12.000 3
## Q2 C07@54.3 C07 54.254 3Fit the model:
## Model-fit
petioleB1.fq <- qtl::fitqtl(Population, pheno.col=7, qtl=twoQTL, get.ests=TRUE, formula=y~Q1+Q2+Q1:Q2, method="hk")## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 180 individuals with missing phenotypes.##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 181
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1 + Q2 + Q1:Q2
##
## df SS MS LOD %var
## Model 8 1147.482 143.43520 9.322536 21.11607
## Error 172 4286.682 24.92257
## Total 180 5434.163
##
##
## Drop one QTL at a time ANOVA table:
## ----------------------------------
## df Type III SS LOD %var F value
## C03@12.0 6 552.5 4.765 10.167 3.69
## C07@54.3 6 648.2 5.535 11.929 4.33
## C03@12.0:C07@54.3 4 198.3 1.777 3.649 1.99
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 29.4103 0.4053 72.556
## C03@12.0a 1.9693 0.6549 3.007
## C03@12.0d 1.1282 0.8365 1.349
## C07@54.3a 2.4536 0.5934 4.135
## C07@54.3d 1.0573 0.8107 1.304
## C03@12.0a:C07@54.3a -1.2842 0.9930 -1.293
## C03@12.0d:C07@54.3a -0.7378 1.2171 -0.606
## C03@12.0a:C07@54.3d 2.1137 1.3099 1.614
## C03@12.0d:C07@54.3d -2.0065 1.6731 -1.199The LOD for interaction is low. The model is reduced to an additive model only. Scan for additional QTL:
## Scan for additional QTL and check
add.qtl <- qtl::addqtl(Population, qtl=twoQTL, method="hk", formula=y~Q1+Q2, pheno.col=7)## Warning in qtl::addqtl(Population, qtl = twoQTL, method = "hk", formula = y ~ : Dropping 180 individuals with missing phenotypes.## chr pos lod pval
## Bn-scaff_16397_1-p165703 C06 47.8 3.53 0.0534Evidence for the QTL on C06 increased from 2.79 to 3.53 (p = 0.053).
Adjust the model and check it:
## Make the model
print(threeQTL <- qtl::makeqtl(Population, chr=c("C03","C07","C06"), pos=c(12,54.2,47.8), what="prob"))## QTL object containing genotype probabilities.
##
## name chr pos n.gen
## Q1 C03@12.0 C03 12.000 3
## Q2 C07@54.3 C07 54.254 3
## Q3 C06@47.8 C06 47.832 3Calculate the penalties:
## main heavy light
## 3.601356 5.866595 3.415472Perform the search:
## Best-model search
print(stepPetioleB1 <- qtl::stepwiseqtl(Population, additive.only=TRUE, qtl=threeQTL, max.qtl=6,
penalties=pen, verbose=TRUE, pheno.col=7, method="hk", keeptrace=TRUE))## -Initial scan
## ---Starting at a model with 3 QTL
## ---Refining positions
## --- Moved a bit
## ** new best ** (pLOD increased by 1.0609)
## no.qtl = 3 pLOD = 1.060851 formula: y ~ Q1+Q2+Q3
## -Step 1
## ---Scanning for additive qtl
## plod = 0.6677621
## ---Refining positions
## no.qtl = 4 pLOD = 0.6677621 formula: y ~ Q1 + Q2 + Q3 + Q4
## -Step 2
## ---Scanning for additive qtl
## plod = -0.3734606
## ---Refining positions
## --- Moved a bit
## no.qtl = 5 pLOD = 0.1394927 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5
## -Step 3
## ---Scanning for additive qtl
## plod = -1.44203
## ---Refining positions
## --- Moved a bit
## no.qtl = 6 pLOD = -1.426128 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q6
## -Starting backward deletion
## ---Dropping Q6
## no.qtl = 5 pLOD = 0.1285551 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5
## ---Refining positions
## --- Moved a bit
## ---Dropping Q5
## no.qtl = 4 pLOD = 0.5683411 formula: y ~ Q1 + Q2 + Q3 + Q4
## ---Refining positions
## --- Moved a bit
## ---Dropping Q4
## no.qtl = 3 pLOD = 1.060851 formula: y ~ Q1 + Q2 + Q3
## ---Refining positions
## ---Dropping Q3
## no.qtl = 2 pLOD = 1.121144 formula: y ~ Q1 + Q2
## ---Refining positions
## ** new best ** (pLOD increased by 0.0603)
## ---Dropping Q1
## no.qtl = 1 pLOD = 0.9563585 formula: y ~ Q1
## ---Refining positions
## ---One last pass through refineqtl
## QTL object containing genotype probabilities.
##
## name chr pos n.gen
## Q1 C03@16.8 C03 16.790 3
## Q2 C07@54.3 C07 54.254 3
##
## Formula: y ~ Q1 + Q2
##
## pLOD: 1.121The best fit is for the QTL on C03 and C07 with additive effects. The position on C03 moved slightly. We adjust the model and check the effect:
## Redefine model
twoQTL <- qtl::makeqtl(Population, chr=c("C03","C07"), pos=c(12, 54.25), what="prob")
reftwoQTL <- qtl::refineqtl(Population, pheno.col=7, method="hk", qtl=twoQTL, keeplodprofile=TRUE)## pos: 12 54.25374
## Iteration 1
## Q2 pos: 54.25374 -> 54.25374
## LOD increase: 0
## Q1 pos: 12 -> 16.79003
## LOD increase: 0.779
## all pos: 12 54.25374 -> 16.79003 54.25374
## LOD increase at this iteration: 0.779
## Iteration 2
## Q1 pos: 16.79003 -> 16.79003
## LOD increase: 0
## Q2 pos: 54.25374 -> 54.25374
## LOD increase: 0
## all pos: 16.79003 54.25374 -> 16.79003 54.25374
## LOD increase at this iteration: 0
## overall pos: 12 54.25374 -> 16.79003 54.25374
## LOD increase overall: 0.779And fit the model:
## Model-fit
petioleB1.final <- qtl::fitqtl(Population, pheno.col=7, qtl=reftwoQTL, get.ests=TRUE, formula=y~Q1+Q2, method="hk")## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 180 individuals with missing phenotypes.##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 181
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1 + Q2
##
## df SS MS LOD %var
## Model 4 1037.164 259.29107 8.323857 19.086
## Error 176 4396.999 24.98295
## Total 180 5434.163
##
##
## Drop one QTL at a time ANOVA table:
## ----------------------------------
## df Type III SS LOD %var F value
## C03@16.8 2 442.2 3.766 8.137 8.85
## C07@54.3 2 542.0 4.568 9.973 10.85
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 29.3915 0.3776 77.832
## C03@16.8a 2.0132 0.5575 3.611
## C03@16.8d 1.2586 0.7580 1.660
## C07@54.3a 2.3568 0.5195 4.536
## C07@54.3d 1.0500 0.7493 1.4014.2.4 Third petiole-assay
Define the model:
## Make the model
print(twoQTL <- qtl::makeqtl(Population, chr=c("C03","C07"), pos=c(12,46.9), what="prob"))## QTL object containing genotype probabilities.
##
## name chr pos n.gen
## Q1 C03@12.0 C03 12.000 3
## Q2 C07@46.9 C07 46.927 3Fit the model:
## Model-fit
petioleB3.fq <- qtl::fitqtl(Population, pheno.col=9, qtl=twoQTL, get.ests=TRUE, formula=y~Q1+Q2, method="hk")## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 192 individuals with missing phenotypes.##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 169
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1 + Q2
##
## df SS MS LOD %var
## Model 4 491.2783 122.8196 4.265716 10.97373
## Error 164 3985.5780 24.3023
## Total 168 4476.8563
##
##
## Drop one QTL at a time ANOVA table:
## ----------------------------------
## df Type III SS LOD %var F value
## C03@12.0 2 168.0 1.515 3.752 3.46
## C07@46.9 2 218.2 1.956 4.874 4.49
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 25.1441 0.3864 65.068
## C03@12.0a 1.5579 0.6128 2.542
## C03@12.0d -0.6346 0.7953 -0.798
## C07@46.9a 1.3006 0.5582 2.330
## C07@46.9d 1.5994 0.7623 2.098## Scan for additional QTL and check
add.qtl <- qtl::addqtl(Population, qtl=twoQTL, method="hk", formula=y~Q1+Q2, pheno.col=9)## Warning in qtl::addqtl(Population, qtl = twoQTL, method = "hk", formula = y ~ : Dropping 192 individuals with missing phenotypes.## There were no LOD peaks above the threshold.Calculate the penalites:
## main heavy light
## 3.563503 6.042247 3.707353Perform the search:
## Best-model search
print(stepPetioleB3 <- qtl::stepwiseqtl(Population, qtl=twoQTL, max.qtl=6,
penalties=pen, verbose=TRUE, pheno.col=9, method="hk",
keeptrace=TRUE))## -Initial scan
## ---Starting at a model with 2 QTL
## ---Refining positions
## --- Moved a bit
## no.qtl = 2 pLOD = -1.35889 formula: y ~ Q1+Q2
## -Step 1
## ---Scanning for additive qtl
## plod = -2.349011
## ---Scanning for QTL interacting with Q1
## plod = -5.329715
## ---Scanning for QTL interacting with Q2
## plod = -4.750527
## ---Look for additional interactions
## plod = -4.320231
## ---Refining positions
## no.qtl = 3 pLOD = -2.349011 formula: y ~ Q1 + Q2 + Q3
## -Step 2
## ---Scanning for additive qtl
## plod = -4.099295
## ---Scanning for QTL interacting with Q1
## plod = -6.349757
## ---Scanning for QTL interacting with Q2
## plod = -5.402249
## ---Scanning for QTL interacting with Q3
## plod = -5.710327
## ---Look for additional interactions
## plod = -5.435507
## ---Refining positions
## no.qtl = 4 pLOD = -4.099295 formula: y ~ Q1 + Q2 + Q3 + Q4
## -Step 3
## ---Scanning for additive qtl
## plod = -5.871496
## ---Scanning for QTL interacting with Q1
## plod = -7.706223
## ---Scanning for QTL interacting with Q2
## plod = -7.490073
## ---Scanning for QTL interacting with Q3
## plod = -6.791686
## ---Scanning for QTL interacting with Q4
## plod = -7.192253
## ---Look for additional interactions
## plod = -5.710327
## ---Refining positions
## --- Moved a bit
## no.qtl = 4 pLOD = -5.667729 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q3:Q4
## -Step 4
## ---Scanning for additive qtl
## plod = -7.038333
## ---Scanning for QTL interacting with Q1
## plod = -8.451984
## ---Scanning for QTL interacting with Q2
## plod = -8.360799
## ---Scanning for QTL interacting with Q3
## plod = -10.70884
## ---Scanning for QTL interacting with Q4
## plod = -11.18362
## ---Look for additional interactions
## plod = -8.750747
## ---Refining positions
## --- Moved a bit
## no.qtl = 5 pLOD = -6.72845 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q3:Q4 + Q5
## -Step 5
## ---Scanning for additive qtl
## plod = -8.093837
## ---Scanning for QTL interacting with Q1
## plod = -10.85961
## ---Scanning for QTL interacting with Q2
## plod = -8.652127
## ---Scanning for QTL interacting with Q3
## plod = -11.87585
## ---Scanning for QTL interacting with Q4
## plod = -12.28263
## ---Scanning for QTL interacting with Q5
## plod = -9.155756
## ---Look for additional interactions
## plod = -10.10689
## ---Refining positions
## --- Moved a bit
## no.qtl = 6 pLOD = -7.935273 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q3:Q4 + Q5 + Q6
## -Starting backward deletion
## ---Dropping Q4
## no.qtl = 5 pLOD = -4.452377 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5
## ---Refining positions
## --- Moved a bit
## ---Dropping Q5
## no.qtl = 4 pLOD = -2.687312 formula: y ~ Q1 + Q2 + Q3 + Q4
## ---Refining positions
## --- Moved a bit
## ---Dropping Q3
## no.qtl = 3 pLOD = -1.43831 formula: y ~ Q1 + Q2 + Q3
## ---Refining positions
## --- Moved a bit
## ---Dropping Q1
## no.qtl = 2 pLOD = -0.5949972 formula: y ~ Q1 + Q2
## ---Refining positions
## --- Moved a bit
## ---Dropping Q2
## no.qtl = 1 pLOD = -1.243681 formula: y ~ Q1
## ---Refining positions
## --- Moved a bit
## Null QTL modelEstimate the effects:
## Fit the model
petioleB3.fq <- qtl::fitqtl(Population, pheno.col=9, qtl=twoQTL, get.ests=TRUE, formula=y~Q1+Q2, method="hk")## Warning in fitqtlengine(pheno = pheno, qtl = qtl, covar = covar, formula = formula, : Dropping 192 individuals with missing phenotypes.##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 169
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1 + Q2
##
## df SS MS LOD %var
## Model 4 491.2783 122.8196 4.265716 10.97373
## Error 164 3985.5780 24.3023
## Total 168 4476.8563
##
##
## Drop one QTL at a time ANOVA table:
## ----------------------------------
## df Type III SS LOD %var F value
## C03@12.0 2 168.0 1.515 3.752 3.46
## C07@46.9 2 218.2 1.956 4.874 4.49
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 25.1441 0.3864 65.068
## C03@12.0a 1.5579 0.6128 2.542
## C03@12.0d -0.6346 0.7953 -0.798
## C07@46.9a 1.3006 0.5582 2.330
## C07@46.9d 1.5994 0.7623 2.0985 QTL-intervals
C01
## chr pos lod
## Bn-scaff_24947_1-p13145 C01 46.24290 4.679284
## Bn-scaff_19564_1-p17934 C01 52.97992 6.229603
## Bn-scaff_18672_1-p178765 C01 57.85777 4.782165## chr pos lod
## Bn-scaff_17827_1-p427361 C01 36.3964 4.744319
## cC01.loc51 C01 51.0000 6.160004
## Bn-scaff_16226_1-p1835 C01 56.4642 4.355473C03
## chr pos lod
## Bn-scaff_16614_1-p2191153 C03 0.00000 2.852941
## cC03.loc12 C03 12.00000 3.787694
## Bn-scaff_15877_1-p293599 C03 18.05588 2.941277## chr pos lod
## Bn-scaff_16614_1-p2191153 C03 0.00000 2.907152
## Bn-scaff_18936_1-p269153 C03 13.12951 3.223499
## Bn-scaff_16135_1-p168458 C03 96.16195 1.537397Note:
Estimated QTL-intervals in the third petiole-assay behave poor due to the low significance of p3QTLb1.
C07
## chr pos lod
## Bn-scaff_16069_1-p3479340 C07 49.15440 3.134250
## Bn-scaff_16069_1-p4306874 C07 54.25374 4.557715
## Bn-scaff_16110_1-p3075966 C07 62.29354 2.760066## chr pos lod
## Bn-scaff_16069_1-p1469430 C07 43.57177 2.444004
## Bn-scaff_16069_1-p2611780 C07 46.92700 3.775532
## Bn-scaff_16110_1-p1330898 C07 75.98592 1.933497## chr pos lod
## Bn-scaff_16110_1-p1938114 C07 71.37864 1.709399
## Bn-scaff_16110_1-p976517 C07 79.55689 4.072145
## Bn-scaff_16110_1-p426547 C07 80.54132 3.694813C06
## chr pos lod
## Bn-scaff_15818_2-p1610442 C06 17.85096 1.275519
## Bn-scaff_16397_1-p165703 C06 47.83219 2.766733
## Bn-scaff_24104_1-p170165 C06 72.20666 1.056143## chr pos lod
## Bn-scaff_15763_1-p1444439 C06 22.90557 0.2574246
## Bn-scaff_16397_1-p360292 C06 46.84870 2.6889376
## Bn-scaff_24104_1-p170165 C06 72.20666 0.25415706 Final plots
6.1 Frequency distribution
Prepare data for plotting:
## Leaf-assay Population A
LeafA <- as.data.frame(Population$pheno[1:180,2])
colnames(LeafA) <- "LeafA"
LeafA <- na.omit(LeafA)
## Plot
DisLeafA <- ggplot2::ggplot(LeafA)+
ggplot2::geom_histogram(aes(LeafA), binwidth=30, color="darkslategray", fill="lightseagreen", na.rm=TRUE)+
ggplot2::geom_vline(aes(xintercept=mean(LeafA)), color="black", linetype="dashed", size=1)+
ggplot2::scale_y_continuous(expand=c(0,0), limits=c(0,25))+
ggplot2::scale_x_continuous(expand=c(0,0), limits=c(0, 1500), breaks=c(0,500,1000,1500))+
ggplot2::labs(title=~1^{st}~"Leaf-assay (Population A)", y = "Frequency", x ="Leaf-lesion area"~"["~mm^{2}~"]")+
ggplot2::theme_light()+
ggplot2::theme(plot.title=element_text(hjust=0.5, size=11),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.ticks.y=element_line(size=0.8, color="black"),
axis.ticks.x=element_line(size=0.8, color="black"),
axis.text.y=element_text(size=11, color="black"),
axis.text.x=element_text(size=11, color="black"),
panel.border=element_rect(linetype="solid", colour="black"),
strip.background=element_rect(colour="black", fill="gray80"),
strip.text.x=element_text(colour="black", face="bold"),
plot.margin = unit(c(0.2,0.4,0.2,0.2), "cm"))
## First leaf-assay Population B
LeafB1 <- as.data.frame(Population$pheno[181:361,4])
colnames(LeafB1) <- "LeafB1"
LeafB1 <- na.omit(LeafB1)
## Plot
DisLeafB1 <- ggplot2::ggplot(LeafB1)+
ggplot2::geom_histogram(aes(LeafB1), binwidth=30, color="darkslategray", fill="lightseagreen", na.rm=TRUE)+
ggplot2::geom_vline(aes(xintercept=mean(LeafB1)), color="black", linetype="dashed", size=1)+
ggplot2::scale_y_continuous(expand=c(0,0), limits=c(0,25))+
ggplot2::scale_x_continuous(expand=c(0,0), limits=c(0, 1500), breaks=c(0,500,1000,1500))+
ggplot2::labs(title=~1^{st}~"Leaf-assay (Population B)", y = "Frequency", x ="Leaf-lesion area"~"["~mm^{2}~"]")+
ggplot2::theme_light()+
ggplot2::theme(plot.title=element_text(hjust=0.5, size=11),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.ticks.y=element_line(size=0.8, color="black"),
axis.ticks.x=element_line(size=0.8, color="black"),
axis.text.y=element_text(size=11, color="black"),
axis.text.x=element_text(size=11, color="black"),
panel.border=element_rect(linetype="solid", colour="black"),
strip.background=element_rect(colour="black", fill="gray80"),
strip.text.x=element_text(colour="black", face="bold"),
plot.margin = unit(c(0.2,0.4,0.2,0.2), "cm"))
## Second leaf-assay Population B
LeafB2 <- as.data.frame(Population$pheno[181:361,5])
colnames(LeafB2) <- "LeafB2"
LeafB2 <- na.omit(LeafB2)
## Plot
DisLeafB2 <- ggplot2::ggplot(LeafB2)+
ggplot2::geom_histogram(aes(LeafB2), binwidth=30, color="darkslategray", fill="lightseagreen", na.rm=TRUE)+
ggplot2::geom_vline(aes(xintercept=mean(LeafB2)), color="black", linetype="dashed", size=1)+
ggplot2::scale_y_continuous(expand=c(0,0), limits=c(0,25))+
ggplot2::scale_x_continuous(expand=c(0,0), limits=c(0, 1500), breaks=c(0,500,1000,1500))+
ggplot2::labs(title=~2^{nd}~"Leaf-assay (Population B)", y = "Frequency", x ="Leaf-lesion area"~"["~mm^{2}~"]")+
ggplot2::theme_light()+
ggplot2::theme(plot.title=element_text(hjust=0.5, size=11),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.ticks.y=element_line(size=0.8, color="black"),
axis.ticks.x=element_line(size=0.8, color="black"),
axis.text.y=element_text(size=11, color="black"),
axis.text.x=element_text(size=11, color="black"),
panel.border=element_rect(linetype="solid", colour="black"),
strip.background=element_rect(colour="black", fill="gray80"),
strip.text.x=element_text(colour="black", face="bold"),
plot.margin = unit(c(0.2,0.4,0.2,0.2), "cm"))
## Third leaf-assay Population B
LeafB3 <- as.data.frame(Population$pheno[181:361,6])
colnames(LeafB3) <- "LeafB3"
LeafB3 <- na.omit(LeafB3)
## Plot
DisLeafB3 <- ggplot2::ggplot(LeafB3)+
ggplot2::geom_histogram(aes(LeafB3), binwidth=30, color="darkslategray", fill="lightseagreen", na.rm=TRUE)+
ggplot2::geom_vline(aes(xintercept=mean(LeafB3)), color="black", linetype="dashed", size=1)+
ggplot2::scale_y_continuous(expand=c(0,0), limits=c(0,25))+
ggplot2::scale_x_continuous(expand=c(0,0), limits=c(0, 1500), breaks=c(0,500,1000,1500))+
ggplot2::labs(title=~3^{rd}~"Leaf-assay (Population B)", y = "Frequency", x ="Leaf-lesion area"~"["~mm^{2}~"]")+
ggplot2::theme_light()+
ggplot2::theme(plot.title=element_text(hjust=0.5, size=11),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.ticks.y=element_line(size=0.8, color="black"),
axis.ticks.x=element_line(size=0.8, color="black"),
axis.text.y=element_text(size=11, color="black"),
axis.text.x=element_text(size=11, color="black"),
panel.border=element_rect(linetype="solid", colour="black"),
strip.background=element_rect(colour="black", fill="gray80"),
strip.text.x=element_text(colour="black", face="bold"),
plot.margin = unit(c(0.2,0.4,0.2,0.2), "cm"))
## Petiole-assay Population A
PetioleA <- as.data.frame(Population$pheno[1:180,3])
colnames(PetioleA) <- "PetioleA"
PetioleA <- na.omit(PetioleA)
## Plot
DisPetioleA <- ggplot2::ggplot(PetioleA, aes(PetioleA))+
ggplot2::geom_histogram(position="identity", binwidth=1, color="darkslategray", fill="lightseagreen", na.rm=TRUE)+
ggplot2::geom_vline(aes(xintercept=mean(PetioleA)), color="black", linetype="dashed", size=1)+
ggplot2::scale_y_continuous(expand=c(0,0), limits=c(0,22))+
ggplot2::scale_x_continuous(expand=c(0,0), limits=c(0,60), breaks=c(10,20,30,40,50,60))+
ggplot2::labs(title=~1^{st}~"Petiole-assay (Population A)", y = "Frequency", x ="Petiole-lesion length [mm]")+
ggplot2::theme_light()+
ggplot2::theme(plot.title=element_text(hjust=0.5, size=11),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.ticks.y=element_line(size=0.8, color="black"),
axis.ticks.x=element_line(size=0.8, color="black"),
axis.text.y=element_text(size=11, color="black"),
axis.text.x=element_text(size=11, color="black"),
panel.border=element_rect(linetype="solid", colour="black"),
strip.background=element_rect(colour="black", fill="gray80"),
strip.text.x=element_text(colour="black", face="bold"),
plot.margin = unit(c(0.2,0.4,0.2,0.2), "cm"))
## First petiole-assay Population B
PetioleB1 <- as.data.frame(Population$pheno[181:361,7])
colnames(PetioleB1) <- "PetioleB1"
PetioleB1 <- na.omit(PetioleB1)
## Plot
DisPetioleB1 <- ggplot2::ggplot(PetioleB1)+
ggplot2::geom_histogram(aes(PetioleB1), binwidth=1, color="darkslategray", fill="lightseagreen", na.rm=TRUE)+
ggplot2::geom_vline(aes(xintercept=mean(PetioleB1)), color="black", linetype="dashed", size=1)+
ggplot2::scale_y_continuous(expand=c(0,0), limits=c(0,22))+
ggplot2::scale_x_continuous(expand=c(0,0), limits=c(0, 60), breaks=c(10,20,30,40,50,60))+
ggplot2::labs(title=~1^{st}~"Petiole-assay (Population B)", y = "Frequency", x ="Petiole-lesion length [mm]")+
ggplot2::theme_light()+
ggplot2::theme(plot.title=element_text(hjust=0.5, size=11),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.ticks.y=element_line(size=0.8, color="black"),
axis.ticks.x=element_line(size=0.8, color="black"),
axis.text.y=element_text(size=11, color="black"),
axis.text.x=element_text(size=11, color="black"),
panel.border=element_rect(linetype="solid", colour="black"),
strip.background=element_rect(colour="black", fill="gray80"),
strip.text.x=element_text(colour="black", face="bold"),
plot.margin = unit(c(0.2,0.4,0.2,0.2), "cm"))
## Second Petiole-assay Population B
PetioleB2 <- as.data.frame(Population$pheno[181:361,8])
colnames(PetioleB2) <- "PetioleB2"
PetioleB2 <- na.omit(PetioleB2)
## Plot
DisPetioleB2 <- ggplot2::ggplot(PetioleB2)+
ggplot2::geom_histogram(aes(PetioleB2), binwidth=1, color="darkslategray", fill="lightseagreen", na.rm=TRUE)+
ggplot2::geom_vline(aes(xintercept=mean(PetioleB2)), color="black", linetype="dashed", size=1)+
ggplot2::scale_y_continuous(expand=c(0,0), limits=c(0,22))+
ggplot2::scale_x_continuous(expand=c(0,0), limits=c(0, 60), breaks=c(10,20,30,40,50,60))+
ggplot2::labs(title=~2^{nd}~"Petiole-assay (Population B)", y = "Frequency", x ="Petiole-lesion length [mm]")+
ggplot2::theme_light()+
ggplot2::theme(plot.title=element_text(hjust=0.5, size=11),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.ticks.y=element_line(size=0.8, color="black"),
axis.ticks.x=element_line(size=0.8, color="black"),
axis.text.y=element_text(size=11, color="black"),
axis.text.x=element_text(size=11, color="black"),
panel.border=element_rect(linetype="solid", colour="black"),
strip.background=element_rect(colour="black", fill="gray80"),
strip.text.x=element_text(colour="black", face="bold"),
plot.margin = unit(c(0.2,0.4,0.2,0.2), "cm"))
## Third Petiole-assay Population B
PetioleB3 <- as.data.frame(Population$pheno[181:361,9])
colnames(PetioleB3) <- "PetioleB3"
PetioleB3 <- na.omit(PetioleB3)
## Plot
DisPetioleB3 <- ggplot2::ggplot(PetioleB3)+
ggplot2::geom_histogram(aes(PetioleB3), binwidth=1, color="darkslategray", fill="lightseagreen", na.rm=TRUE)+
ggplot2::geom_vline(aes(xintercept=mean(PetioleB3)), color="black", linetype="dashed", size=1)+
ggplot2::scale_y_continuous(expand=c(0,0), limits=c(0,22))+
ggplot2::scale_x_continuous(expand=c(0,0), limits=c(0, 60), breaks=c(10,20,30,40,50,60))+
ggplot2::labs(title=~3^{rd}~"Petiole-assay (Population B)", y = "Frequency", x ="Petiole-lesion length [mm]")+
ggplot2::theme_light()+
ggplot2::theme(plot.title=element_text(hjust=0.5, size=11),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.ticks.y=element_line(size=0.8, color="black"),
axis.ticks.x=element_line(size=0.8, color="black"),
axis.text.y=element_text(size=11, color="black"),
axis.text.x=element_text(size=11, color="black"),
panel.border=element_rect(linetype="solid", colour="black"),
strip.background=element_rect(colour="black", fill="gray80"),
strip.text.x=element_text(colour="black", face="bold"),
plot.margin = unit(c(0.2,0.4,0.2,0.2), "cm"))Combine plots:
## Bring all together
PlotDistributions <- cowplot::plot_grid(DisLeafA, DisLeafB1, DisLeafB2, DisLeafB3,
DisPetioleA, DisPetioleB1, DisPetioleB2, DisPetioleB3,
labels=c("A","B","C","D","E","F","G","H"),
ncol=4, nrow=2, align = "hv", label_size = 16)
dev.new(width=20, height=18)Plot distributions:
Frequency distributions
## For export
#png("./pyQTL/Abbildungen/Verteilungen.png", res=280, width=3000, height=2000, type="cairo")
#PlotDistributions
#dev.off()6.2 LOD-profiles
Prepare data for plotting:
## Convert LOD-profile to data frame
## Petioles
LODPetioleA <- as.data.frame(out.petioleA.hk)
LODPetioleB1 <- as.data.frame(out.petioleB1.hk)
LODPetioleB2 <- as.data.frame(out.petioleB2.hk)
LODPetioleB3 <- as.data.frame(out.petioleB3.np)
## Peaks
maxLODPetioleA <- as.data.frame(max(out.petioleA.hk))
max1LODPetioleB1 <- as.data.frame(max(out.petioleB1.hk[out.petioleB1.hk$chr=="C03",]))
max2LODPetioleB1 <- as.data.frame(max(out.petioleB1.hk[out.petioleB1.hk$chr=="C07",]))
max1LODPetioleB3 <- as.data.frame(max(out.petioleB3.np[out.petioleB3.np$chr=="C03",]))
max2LODPetioleB3 <- as.data.frame(max(out.petioleB3.np[out.petioleB3.np$chr=="C07",]))
## Plot for single-scan
labs <- c('C01' = "Chromosome 1",
'C02' = "Chromosome 2",
'C03' = "Chromosome 3",
'C04a' = "Chromosome 4a",
'C04b' = "Chromosome 4b",
'C05' = "Chromosome 5",
'C06' = "Chromosome 6",
'C07' = "Chromosome 7",
'C08' = "Chromosome 8",
'C09' = "Chromosome 9")
##
LODPetioles <- ggplot2::ggplot(LODPetioleA) +
ggplot2::geom_line(aes(x=pos, y=lod), size=0.5, alpha=0.8, color="red")+
ggplot2::geom_line(aes(x=pos, y=LODPetioleB1$lod), size=0.5, alpha=0.8, color="darkslategray")+
ggplot2::geom_line(aes(x=pos, y=LODPetioleB2$lod), size=0.5, alpha=0.8, color="deepskyblue")+
ggplot2::geom_line(aes(x=pos, y=LODPetioleB3$lod), size=0.5, alpha=0.8, color="lightseagreen")+
ggplot2::geom_point(data = maxLODPetioleA, aes(x=pos, y=lod))+
ggplot2::geom_point(data = max1LODPetioleB1, aes(x=pos, y=lod))+
ggplot2::geom_point(data = max2LODPetioleB1, aes(x=pos, y=lod))+
ggplot2::geom_point(data = max1LODPetioleB3, aes(x=pos, y=lod))+
ggplot2::geom_point(data = max2LODPetioleB3, aes(x=pos, y=lod))+
ggplot2::labs(title="", y="LOD", x="Interval position [cM]")+
ggplot2::scale_y_continuous(expand=c(0,0), breaks=c(0,2,4,6), limits=c(0,7))+
ggplot2::geom_abline(aes(intercept=3.65, slope=0), linetype="solid", size=0.2, colour="red")+
ggplot2::geom_abline(aes(intercept=3.56, slope=0), linetype="solid", size=0.2, colour="darkslategray")+
ggplot2::geom_abline(aes(intercept=3.66, slope=0), linetype="solid", size=0.2, colour="deepskyblue")+
ggplot2::geom_abline(aes(intercept=3.49, slope=0), linetype="solid", size=0.2, colour="lightseagreen")+
ggrepel::geom_text_repel(data = maxLODPetioleA, aes(x = pos, y = lod, label = "pQTLA"), size=2.5, nudge_x = -25, nudge_y = 0.5)+
ggrepel::geom_text_repel(data = max1LODPetioleB1, aes(x = pos, y = lod, label = "p1QTLb1"), size=2.5, nudge_x = -35, nudge_y = 0.7)+
ggrepel::geom_text_repel(data = max2LODPetioleB1, aes(x = pos, y = lod, label = "p1QTLb2"), size=2.5, nudge_x = 25, nudge_y = 0.7)+
ggrepel::geom_text_repel(data = max1LODPetioleB3, aes(x = pos, y = lod, label = "p3QTLb1"), size=2.5, nudge_x = 55, nudge_y = 0.1)+
ggrepel::geom_text_repel(data = max2LODPetioleB3, aes(x = pos, y = lod, label = "p3QTLb2"), size=2.5, nudge_x = -25, nudge_y = 0.4)+
ggplot2::facet_wrap(~ chr, scales = "free_x", ncol=10, labeller=as_labeller(labs))+
ggplot2::theme_light()+
ggplot2::theme(plot.title=element_text(hjust=0.5, size=14),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.ticks.y=element_line(size=0.8, color="black"),
axis.ticks.x=element_line(size=0.8, color="black"),
axis.text.y=element_text(size=11, color="black"),
axis.text.x=element_text(size=11, color="black"),
panel.border=element_rect(linetype="solid", colour="black"),
strip.background=element_rect(colour="black", fill="gray80"),
strip.text.x=element_text(colour="black", face="bold"))
## Leaves
LODLeafA <- as.data.frame(out.leafA.hk)
LODLeafB1 <- as.data.frame(out.leafB1.hk)
LODLeafB2 <- as.data.frame(out.leafB2.hk)
LODLeafB3 <- as.data.frame(out.leafB3.hk)
## Peaks
maxLODLeafB2 <- as.data.frame(max(out.leafB2.hk[out.leafB2.hk$chr=="C01",]))
maxLODLeafB3 <- as.data.frame(max(out.leafB3.hk[out.leafB3.hk$chr=="C07",]))
##
LODLeaves <- ggplot2::ggplot(LODLeafA) +
ggplot2::geom_line(aes(x=pos, y=lod), size=0.5, alpha=0.8, color="red")+
ggplot2::geom_line(aes(x=pos, y=LODLeafB1$lod), size=0.5, alpha=0.8, color="darkslategray")+
ggplot2::geom_line(aes(x=pos, y=LODLeafB2$lod), size=0.5, alpha=0.8, color="deepskyblue")+
ggplot2::geom_line(aes(x=pos, y=LODLeafB3$lod), size=0.5, alpha=0.8, color="lightseagreen")+
ggplot2::geom_point(data = maxLODLeafB2, aes(x=pos, y=lod))+
ggplot2::geom_point(data = maxLODLeafB3, aes(x=pos, y=lod))+
ggplot2::labs(title="", y="LOD", x="Interval position [cM]")+
ggplot2::scale_y_continuous(expand=c(0,0), breaks=c(0,2,4,6), limits=c(0,7))+
ggplot2::geom_abline(aes(intercept=3.61, slope=0), linetype="solid", size=0.2, colour="red")+
ggplot2::geom_abline(aes(intercept=3.63, slope=0), linetype="solid", size=0.2, colour="darkslategray")+
ggplot2::geom_abline(aes(intercept=3.64, slope=0), linetype="solid", size=0.2, colour="deepskyblue")+
ggplot2::geom_abline(aes(intercept=3.59, slope=0), linetype="solid", size=0.2, colour="lightseagreen")+
ggrepel::geom_text_repel(data = maxLODLeafB2, aes(x = pos, y = lod, label = "l2QTLb"), size=2.5, nudge_x = -20, nudge_y = 0.4)+
ggrepel::geom_text_repel(data = maxLODLeafB3, aes(x = pos, y = lod, label = "l3QTLb"), size=2.5, nudge_x = -25, nudge_y = 0.4)+
ggplot2::facet_wrap(~ chr, scales = "free_x", ncol=10, labeller=as_labeller(labs))+
ggplot2::theme_light()+
ggplot2::theme(plot.title=element_text(hjust=0.5, size=14),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.ticks.y=element_line(size=0.8, color="black"),
axis.ticks.x=element_line(size=0.8, color="black"),
axis.text.y=element_text(size=11, color="black"),
axis.text.x=element_text(size=11, color="black"),
panel.border=element_rect(linetype="solid", colour="black"),
strip.background=element_rect(colour="black", fill="gray80"),
strip.text.x=element_text(colour="black", face="bold"))Combine plots:
## Bring all together
PlotLODs <- cowplot::plot_grid(LODLeaves, LODPetioles, labels=c("A","B"),
ncol=1, nrow=2, align = "hv", label_size = 20)
dev.new(width=12, height=8)LOD-profiles:
LOD-profiles for Sclerotini-resistance. Top: Leaf-assays. Bottom: Petiole-assays. Red: Scans in Population A; blue to turquoise: scans in Population B
## For export
#png("./pyQTL/Abbildungen/LOD-Profiles.png", res=250, width=3000, height=2000, type="cairo")
#PlotLODs
#dev.off()Finished
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=German_Germany.1252 LC_CTYPE=German_Germany.1252
## [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C
## [5] LC_TIME=German_Germany.1252
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] multcomp_1.4-17 TH.data_1.0-10 MASS_7.3-51.5
## [4] survival_3.1-8 mvtnorm_1.1-1 gplots_3.1.1
## [7] gdata_2.18.0 ggrepel_0.9.1 qtl_1.48-1
## [10] extrafont_0.17 VennDiagram_1.6.20 futile.logger_1.4.3
## [13] ggpubr_0.4.0 lsmeans_2.30-0 emmeans_1.6.0
## [16] gridExtra_2.3 wesanderson_0.3.6 RColorBrewer_1.1-2
## [19] ggsignif_0.6.1 dplyr_1.0.5 cowplot_1.1.1
## [22] ggplot2_3.3.3 readxl_1.3.1
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-144 bitops_1.0-7 tools_3.6.3
## [4] backports_1.2.1 DT_0.18 utf8_1.1.4
## [7] R6_2.5.0 KernSmooth_2.23-16 DBI_1.1.1
## [10] mgcv_1.8-31 colorspace_2.0-0 withr_2.4.2
## [13] tidyselect_1.1.1 curl_4.3 compiler_3.6.3
## [16] extrafontdb_1.0 formatR_1.9 sandwich_3.0-0
## [19] labeling_0.4.2 caTools_1.18.2 scales_1.1.1
## [22] stringr_1.4.0 digest_0.6.27 foreign_0.8-75
## [25] rmarkdown_2.7 rio_0.5.26 pkgconfig_2.0.3
## [28] htmltools_0.5.1.1 highr_0.9 htmlwidgets_1.5.3
## [31] rlang_0.4.10 generics_0.1.0 farver_2.1.0
## [34] jsonlite_1.7.2 zoo_1.8-9 crosstalk_1.1.1
## [37] gtools_3.8.2 zip_2.1.1 car_3.0-10
## [40] magrittr_2.0.1 Matrix_1.2-18 Rcpp_1.0.6
## [43] munsell_0.5.0 fansi_0.4.2 abind_1.4-5
## [46] lifecycle_1.0.0 stringi_1.5.3 yaml_2.2.1
## [49] carData_3.0-4 plyr_1.8.6 parallel_3.6.3
## [52] forcats_0.5.1 crayon_1.4.1 lattice_0.20-38
## [55] haven_2.4.1 splines_3.6.3 hms_1.0.0
## [58] knitr_1.33 pillar_1.6.0 estimability_1.3
## [61] codetools_0.2-16 futile.options_1.0.1 glue_1.4.2
## [64] evaluate_0.14 lambda.r_1.2.4 data.table_1.14.0
## [67] vctrs_0.3.6 Rttf2pt1_1.3.8 cellranger_1.1.0
## [70] gtable_0.3.0 purrr_0.3.4 tidyr_1.1.3
## [73] assertthat_0.2.1 xfun_0.21 openxlsx_4.2.3
## [76] xtable_1.8-4 broom_0.7.6 rstatix_0.7.0
## [79] tibble_3.0.4 ellipsis_0.3.1