QTL - Composite Interval Mapping and (IM)
P. A. Velasquez Vasconez, Mayara Salvian, Vamilton Franzo
16 de junio de 2017
QTL Mapping in mimulus by Composite Interval Mapping and Interval Mapping
We need to export the list of maps using the function write_map let us export the 14 linkage groups previously constructed (https://rpubs.com/PAVelasquezVasconez/280676), and to a file named “maps_list.map”
library("qtl")
library("onemap")
raw_file <- paste(system.file("extdata", package = "onemap"),
"m_feb06.raw", sep = "\\")
fake_f2_qtl <- read.cross("mm", file = "C:\\Users\\alex\\Desktop\\bio\\m_feb06.raw",
mapfile = "maps_list.map")## --Read the following data:
## Type of cross: f2
## Number of individuals: 287
## Number of markers: 418
## Number of phenotypes: 16
## --Cross type: f2We can to comparate our map with output of R/qtl
fake_f2_qtl <- read.cross("mm", file = "C:\\Users\\alex\\Desktop\\bio\\m_feb06.raw",
mapfile = "maps_list.map")## --Read the following data:
## Type of cross: f2
## Number of individuals: 287
## Number of markers: 418
## Number of phenotypes: 16
## --Cross type: f2fake_f2_qtl## 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: 360
## No. markers: 18 22 22 35 19 45 14 32 20 34 21 24 21 33
## Percent genotyped: 89.5
## Genotypes (%): AA:16.8 AB:25.7 BB:22.5 not BB:18.3 not AA:16.8newmap <- est.map(fake_f2_qtl, tol = 1e-6, map.function = "kosambi")
plot.map(fake_f2_qtl, newmap)
We can see minor differences in the distances between markers, in the chromosomes. We can see our map with 14 chromosomes
plot(newmap,xlab="Chromosome",
main="Genetic Map in Mimulus species")
####We investigate the genetic basis of the anther sterility
plot.pheno<-plot.pheno(fake_f2_qtl, pheno.col = 8)
plot(plot.pheno, col="red", xlab="Pollen viability", main="Pollen viability")
fake_f2_qtl <- calc.genoprob(fake_f2_qtl, step = 2)
out.em <- scanone(fake_f2_qtl, pheno.col=8, method = "em")## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 12 individuals with missing phenotypes.out.hk <- scanone(fake_f2_qtl, pheno.col =8, method = "hk")## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 12 individuals with missing phenotypes.plot(out.em, out.hk, col = c("blue", "gray"), lwd=1)
####The infortations for the two models
summary(out.em)## chr pos lod
## c1.loc176 1 176.0 1.959
## c2.loc26 2 26.0 1.804
## c3.loc52 3 52.0 0.386
## c4.loc140 4 140.0 1.445
## c5.loc40 5 40.0 1.364
## CA283 6 162.7 4.064
## c7.loc22 7 22.0 1.484
## MgSTS31 8 91.2 2.104
## c9.loc64 9 64.0 1.715
## c10.loc26 10 26.0 1.760
## c11.loc56 11 56.0 2.226
## BA175 12 136.2 0.817
## c13.loc58 13 58.0 9.150
## c14.loc208 14 208.0 4.152summary(out.hk)## chr pos lod
## c1.loc176 1 176.0 1.964
## c2.loc50 2 50.0 1.808
## c3.loc52 3 52.0 0.379
## c4.loc140 4 140.0 1.463
## c5.loc40 5 40.0 1.371
## c6.loc162 6 162.0 4.130
## c7.loc22 7 22.0 1.491
## MgSTS31 8 91.2 2.114
## c9.loc64 9 64.0 1.723
## c10.loc28 10 28.0 1.757
## c11.loc54 11 54.0 2.057
## BA175 12 136.2 0.822
## c13.loc58 13 58.0 8.311
## c14.loc210 14 210.0 4.052operm.hk <- scanone(fake_f2_qtl, 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 ...MI<-summary(operm.hk)[1,1]
plor.hk<-plot(operm.hk, col=c("blue"))
abline(plor.hk, v=operm.hk[1,1],operm.hk, col="red", lty=6)
plot(out.em, out.hk, col = c("blue", "gray"), lwd=1)
abline(h=operm.hk[1,1], col="red", lty=6)
(QTL_IM<-summary(out.hk, perms=operm.hk, alpha=0.05, pvalues=T))## chr pos lod pval
## c6.loc162 6 162 4.13 0.030
## c13.loc58 13 58 8.31 0.000
## c14.loc210 14 210 4.05 0.036Localization do QTL chromossomes
(lodint(out.hk, chr=6, expandtomarkers=T))## chr pos lod
## MgSTS220 6 125.8921 2.360086
## c6.loc162 6 162.0000 4.129576
## CC126 6 176.6469 1.997829(lodint(out.hk, chr=13, expandtomarkers=T))## chr pos lod
## MgSTS104 13 42.28920 5.844983
## c13.loc58 13 58.00000 8.311105
## CA315 13 81.83048 2.945060(lodint(out.hk, chr=14, expandtomarkers=T))## chr pos lod
## MgSTS494 14 189.9540 1.795480
## c14.loc210 14 210.0000 4.051774
## BC542 14 255.9996 1.685087qc <- c("6", "13", "14")
qp <- c(162, 58, 210)
fake.f2 <- subset(fake_f2_qtl, chr=qc)
fake.f2 <- sim.geno(fake.f2, n.draws=8, step=2, err=0.001)
(qtl <- makeqtl(fake.f2, qc, qp, what="draws"))## QTL object containing imputed genotypes, with 8 imputations.
##
## name chr pos n.gen
## Q1 6@162.0 6 162 3
## Q2 13@58.0 13 58 3
## Q3 14@210.0 14 210 3plot(qtl)
##Interval Mapping
GL<-summary(out.hk, perms=operm.hk, alpha=0.05, pvalues=T)[,1]
Position<-summary(out.hk, perms=operm.hk, alpha=0.05, pvalues=T)[,2]
LOD<-summary(out.hk, perms=operm.hk, alpha=0.05, pvalues=T)[,3]
p.value<-summary(out.hk, perms=operm.hk, alpha=0.05, pvalues=T)[,4]Markers with QTLs
(markers <- find.marker(fake_f2_qtl, c(6,13,14), c(162,58,210)))## [1] "CA283" "CB55" "MgSTS16"geno <- pull.geno(fake_f2_qtl)[,markers]
phe<-find.pheno(fake_f2_qtl,"pv")
phenotipo <- pull.pheno(fake_f2_qtl)[,phe]
fenmarc<-data.frame(cbind(phenotipo,geno))
Effect<-(lm(phenotipo~geno,data=fenmarc))$coefficients[2:4]
(R_1<-summary(lm(phenotipo~geno[,1],data=fenmarc))$r.square)## [1] 0.01212717(R_2<-summary(lm(phenotipo~geno[,2],data=fenmarc))$r.square)## [1] 0.03361422(R_3<-summary(lm(phenotipo~geno[,3],data=fenmarc))$r.square)## [1] 0.0461275R2<-c(R_1,R_2,R_3)Table 2. Linkage group, positions, LOD, effect and R2
(Table<-data.frame("Linkage Group"= c(QTL_IM$chr[1],
QTL_IM$chr[2], QTL_IM$chr[3]),
"Position"=c(QTL_IM$pos[1],
QTL_IM$pos[2], QTL_IM$pos[3]), "LOD"= c(QTL_IM$lod[1],
QTL_IM$lod[2], QTL_IM$lod[3]), "Effect"=Effect,
"R2"= R2))## Linkage.Group Position LOD Effect R2
## genoCA283 6 162 4.129576 0.03369169 0.01212717
## genoCB55 13 58 8.311105 -0.06110671 0.03361422
## genoMgSTS16 14 210 4.051774 0.04990158 0.04612750Composite Interval Mapping (CIM)
qtl_prob <- calc.genoprob(fake_f2_qtl, step=1,
error.prob=0.001, map.function=c("kosambi"),
stepwidth=c("fixed"))
qtl_cim<- sim.geno(qtl_prob, n.draws=16, step=1,
off.end=0, error.prob=0.001,
map.function=c("kosambi"),
stepwidth=c("fixed"))
cov<-(stepwiseqtl(qtl_prob, pheno.col=8, method=c("hk"),
covar = NULL, model=c("normal"),
incl.markers=TRUE, refine.locations=TRUE,
additive.only=FALSE, scan.pairs=FALSE,
keeplodprofile=TRUE, keeptrace=FALSE,
verbose=TRUE, require.fullrank=FALSE))## -Initial scan
## initial lod: 8.279218
## ** new best ** (pLOD increased by 4.7592)
## no.qtl = 1 pLOD = 4.759218 formula: y ~ Q1
## -Step 1
## ---Scanning for additive qtl
## plod = 7.53468
## ---Scanning for QTL interacting with Q1
## plod = 10.15381
## ---Refining positions
## --- Moved a bit
## no.qtl = 2 pLOD = 10.86002 formula: y ~ Q1 + Q2 + Q1:Q2
## ** new best ** (pLOD increased by 6.1008)
## -Step 2
## ---Scanning for additive qtl
## plod = 10.91642
## ---Scanning for QTL interacting with Q1
## plod = 7.826683
## ---Scanning for QTL interacting with Q2
## plod = 7.992514
## ---Look for additional interactions
## ---Refining positions
## --- Moved a bit
## no.qtl = 3 pLOD = 11.11771 formula: y ~ Q1 + Q2 + Q1:Q2 + Q3
## ** new best ** (pLOD increased by 0.2577)
## -Step 3
## ---Scanning for additive qtl
## plod = 11.06593
## ---Scanning for QTL interacting with Q1
## plod = 8.34772
## ---Scanning for QTL interacting with Q2
## plod = 8.803038
## ---Scanning for QTL interacting with Q3
## plod = 9.714898
## ---Look for additional interactions
## plod = 7.816965
## ---Refining positions
## --- Moved a bit
## no.qtl = 4 pLOD = 11.08155 formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4
## -Step 4
## ---Scanning for additive qtl
## plod = 9.944729
## ---Scanning for QTL interacting with Q1
## plod = 8.59307
## ---Scanning for QTL interacting with Q2
## plod = 7.511044
## ---Scanning for QTL interacting with Q3
## plod = 9.163309
## ---Scanning for QTL interacting with Q4
## plod = 9.307072
## ---Look for additional interactions
## plod = 8.805918
## ---Refining positions
## --- Moved a bit
## no.qtl = 5 pLOD = 10.34722 formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4 + Q5
## -Step 5
## ---Scanning for additive qtl
## plod = 9.642601
## ---Scanning for QTL interacting with Q1
## plod = 7.570209
## ---Scanning for QTL interacting with Q2
## plod = 6.837968
## ---Scanning for QTL interacting with Q3
## plod = 9.110951
## ---Scanning for QTL interacting with Q4
## plod = 8.663964
## ---Scanning for QTL interacting with Q5
## plod = 10.07941
## ---Look for additional interactions
## plod = 8.539705
## ---Refining positions
## --- Moved a bit
## no.qtl = 6 pLOD = 10.14991 formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4 + Q5 + Q6 + Q5:Q6
## -Step 6
## ---Scanning for additive qtl
## plod = 9.489171
## ---Scanning for QTL interacting with Q1
## plod = 8.131762
## ---Scanning for QTL interacting with Q2
## plod = 6.871478
## ---Scanning for QTL interacting with Q3
## plod = 8.185784
## ---Scanning for QTL interacting with Q4
## plod = 8.29245
## ---Scanning for QTL interacting with Q5
## plod = 7.015311
## ---Scanning for QTL interacting with Q6
## plod = 7.367391
## ---Look for additional interactions
## plod = 8.683088
## ---Refining positions
## --- Moved a bit
## no.qtl = 7 pLOD = 9.58532 formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4 + Q5 + Q6 + Q5:Q6 + Q7
## -Step 7
## ---Scanning for additive qtl
## plod = 8.730081
## ---Scanning for QTL interacting with Q1
## plod = 5.105551
## ---Scanning for QTL interacting with Q2
## plod = 7.231426
## ---Scanning for QTL interacting with Q3
## plod = 8.188776
## ---Scanning for QTL interacting with Q4
## plod = 7.933156
## ---Scanning for QTL interacting with Q5
## plod = 6.469789
## ---Scanning for QTL interacting with Q6
## plod = 6.297228
## ---Scanning for QTL interacting with Q7
## plod = 9.137022
## ---Look for additional interactions
## plod = 8.303687
## ---Refining positions
## --- Moved a bit
## no.qtl = 8 pLOD = 9.500639 formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4 + Q5 + Q6 + Q5:Q6 + Q7 + Q8 + Q7:Q8
## -Step 8
## ---Scanning for additive qtl
## plod = 7.963167
## ---Scanning for QTL interacting with Q1
## plod = 5.297642
## ---Scanning for QTL interacting with Q2
## plod = 6.030815
## ---Scanning for QTL interacting with Q3
## plod = 8.752074
## ---Scanning for QTL interacting with Q4
## plod = 8.067517
## ---Scanning for QTL interacting with Q5
## plod = 5.702456
## ---Scanning for QTL interacting with Q6
## plod = 6.24545
## ---Scanning for QTL interacting with Q7
## plod = 5.113021
## ---Scanning for QTL interacting with Q8
## plod = 6.500578
## ---Look for additional interactions
## plod = 7.675149
## ---Refining positions
## --- Moved a bit
## no.qtl = 9 pLOD = 8.940112 formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4 + Q5 + Q6 + Q5:Q6 + Q7 + Q8 + Q7:Q8 + Q9 + Q3:Q9
## -Step 9
## ---Scanning for additive qtl
## plod = 8.446806
## ---Scanning for QTL interacting with Q1
## plod = 5.167428
## ---Scanning for QTL interacting with Q2
## plod = 5.039761
## ---Scanning for QTL interacting with Q3
## plod = 6.046117
## ---Scanning for QTL interacting with Q4
## plod = 7.39289
## ---Scanning for QTL interacting with Q5
## plod = 5.719645
## ---Scanning for QTL interacting with Q6
## plod = 5.205946
## ---Scanning for QTL interacting with Q7
## plod = 5.036799
## ---Scanning for QTL interacting with Q8
## plod = 5.540829
## ---Scanning for QTL interacting with Q9
## plod = 5.647068
## ---Look for additional interactions
## plod = 5.716087
## ---Refining positions
## --- Moved a bit
## no.qtl = 10 pLOD = 8.806342 formula: y ~ Q1 + Q2 + Q1:Q2 + Q3 + Q4 + Q5 + Q6 + Q5:Q6 + Q7 + Q8 + Q7:Q8 + Q9 + Q3:Q9 + Q10
## -Starting backward deletion
## ---Dropping Q10
## no.qtl = 9 pLOD = 8.712053 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q7 + Q8 + Q9 + Q1:Q2 + Q5:Q6 + Q7:Q8 + Q3:Q9
## ---Refining positions
## --- Moved a bit
## ---Dropping Q9
## no.qtl = 8 pLOD = 9.15033 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q7 + Q8 + Q1:Q2 + Q5:Q6 + Q7:Q8
## ---Refining positions
## --- Moved a bit
## ---Dropping Q8
## no.qtl = 7 pLOD = 9.186153 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q7 + Q1:Q2 + Q5:Q6
## ---Refining positions
## --- Moved a bit
## ---Dropping Q7
## no.qtl = 6 pLOD = 10.04729 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q1:Q2 + Q5:Q6
## ---Refining positions
## --- Moved a bit
## ---Dropping Q5
## no.qtl = 5 pLOD = 10.20379 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q5 + Q1:Q2
## ---Refining positions
## --- Moved a bit
## ---Dropping Q5
## no.qtl = 4 pLOD = 11.02864 formula: y ~ Q1 + Q2 + Q3 + Q4 + Q1:Q2
## ---Refining positions
## --- Moved a bit
## ---Dropping Q4
## no.qtl = 3 pLOD = 10.98898 formula: y ~ Q1 + Q2 + Q3 + Q1:Q2
## ---Refining positions
## --- Moved a bit
## ---Dropping Q3
## no.qtl = 2 pLOD = 10.51275 formula: y ~ Q1 + Q2 + Q1:Q2
## ---Refining positions
## --- Moved a bit
## ---Dropping Q1:Q2
## no.qtl = 2 pLOD = 6.365183 formula: y ~ Q1 + Q2
## ---Refining positions
## --- Moved a bit
## ---Dropping Q2
## no.qtl = 1 pLOD = 4.759218 formula: y ~ Q1
## ---Refining positions
## ---One last pass through refineqtlplot(cov)
out.cim <- cim(qtl_prob, pheno.col=8, n.marcovar= 4,
window=15, method="hk", error.prob=0.001,
map.function=c("kosambi"))
CIMthresold <- cim(qtl_cim,pheno.col=8,n.marcovar=4,method="hk",
imp.method="imp", error.prob=0.0001,
map.function="kosambi",n.perm=1000)## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : addcovar appears to be over-specified; consider dropping columns.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : addcovar appears to be over-specified; consider dropping columns.
## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : addcovar appears to be over-specified; consider dropping columns.CIM<-summary(CIMthresold)[1,1]
plot_hk_CIM<-plot(CIMthresold, col=c("blue"))
abline(plot_hk_CIM, v=CIMthresold[1,1],CIMthresold, col="red", lty=6)
plot(out.cim, out.hk, out.em,
col = c("blue", "gray", " dark red" ), lwd=0.1)
add.cim.covar(out.cim)
abline(h=CIMthresold[1,1], col="blue", lty=6, )
abline(h=operm.hk[1,1], col="red", lty=6)
####QTL composite Interval Mapping chromosome 6
plot(out.cim,chr=6, out.hk, out.em,
main=substitute("QTL detected in chromossome 6 by CIM and IM"),
col = c("blue", "red", "black"), lwd=0.5)
abline(h=CIMthresold[1,1], col="blue", lty=6)
abline(h=operm.hk[1,1], col="red", lty=6)
####QTL Composite Interval Mapping chromossome 8
plot(out.cim,chr=c(7), out.hk,
main=substitute("QTL detected in chromossome 8 only by CIM"),
col = c("blue", "black"),
xlab="Chromossome 8",lwd=1)
abline(h=CIMthresold[1,1], col="blue", lty=6)
abline(h=operm.hk[1,1], col="black", lty=6)
####QTL Composite Interval Mapping chromossome 13
plot(out.cim, chr=c(13), out.hk,
main=substitute("QTL detected in chromossome 13 by CIM and IM"),
col = c("blue", "black"), lwd=1)
abline(h=CIMthresold[1,1], col="blue", lty=6)
abline(h=operm.hk[1,1], col="black", lty=6)
####QTL Composite Interval Mapping chromossome 14
plot(out.cim,chr=c(14), out.hk,
main=substitute("QTL detected in chromossome 14 by CIM and IM"),
col = c("blue", "gray"), lwd=1)
abline(h=CIMthresold[1,1], col="blue", lty=6)
abline(h=operm.hk[1,1], col="red", lty=6)
####QTLs detection in our genetic map
chr1<-(lodint(out.cim, chr=6, expandtomarkers=T))
chr2<-(lodint(out.cim, chr=8, expandtomarkers=T))
chr3<-(lodint(out.cim, chr=13, expandtomarkers=T))
chr4<-(lodint(out.cim, chr=14, expandtomarkers=T))
qc <- c("6", "8", "13", "14")
qp <- c(135, 91.18, 57, 208)
fake.f2_CIM <- subset(fake_f2_qtl, chr=qc)
fake.f2 <- sim.geno(fake.f2_CIM, n.draws=8, step=2, err=0.001)
(qtl_CIM <- makeqtl(fake.f2, qc, qp, what="draws"))## QTL object containing imputed genotypes, with 8 imputations.
##
## name chr pos n.gen
## Q1 6@134.0 6 134.000 3
## Q2 8@91.2 8 91.177 3
## Q3 13@56.0 13 56.000 3
## Q4 14@208.0 14 208.000 3plot(qtl_CIM)
(markers <- find.marker(fake_f2_qtl, c("6", "8", "13", "14"), c(135, 91.18, 57, 208)))## [1] "BD169" "MgSTS31" "CB55" "MgSTS16"geno <- pull.geno(fake_f2_qtl)[,markers]
phe<-find.pheno(fake_f2_qtl,"pv")
phenotipo <- pull.pheno(fake_f2_qtl)[,phe]
fenmarc2<-data.frame(cbind(phenotipo,geno))
Effect_CIM<-(lm(phenotipo~geno,data=fenmarc2))$coefficients[2:5]
(R_1<-summary(lm(phenotipo~geno[,1],data=fenmarc2))$r.square)## [1] 0.04848043(R_2<-summary(lm(phenotipo~geno[,2],data=fenmarc2))$r.square)## [1] 0.01941236(R_3<-summary(lm(phenotipo~geno[,3],data=fenmarc2))$r.square)## [1] 0.03361422(R_4<-summary(lm(phenotipo~geno[,4],data=fenmarc2))$r.square)## [1] 0.0461275R.2<-c(R_1,R_2,R_3,R_4)Table 2. Linkage group, positions, LOD, effect and R2
(Table<-data.frame("Linkage Group"= c(chr1[1,1],chr2[1,1],chr3[1,1],
chr4[1,1]),
"Position"=c(chr1[2,2],chr2[2,2],chr3[2,2],
chr4[2,2]),
"LOD"= c(chr1[2,3],chr2[2,3],chr3[2,3],
chr4[2,3]),
"Effect"=Effect_CIM,
"R2"= R.2))## Linkage.Group Position LOD Effect R2
## genoBD169 6 160 6.293950 -0.03157624 0.04848043
## genoMgSTS31 8 148 2.682043 0.04344111 0.01941236
## genoCB55 13 58 10.898914 -0.08172548 0.03361422
## genoMgSTS16 14 216 3.014661 0.03251362 0.04612750