Lovell — Jul 14, 2014, 6:33 PM
#############################
### Title: Re-analysis of Drydown QTLs
### Author: JTLovell
### Date Created: 8-August 2013
### Modifications (Date):
#############################
###Notes:###
#############################
# In the following code, I use three main objects: "cross","grid","mod"
# these are the input data, pseudomarker grid and the qtl models
#############################
###Libraries:###
#############################
library(qtl)
library(snow)
library(rlecuyer)
library(plyr)
library(ggplot2)
options(warn=-1)
sessionInfo()
R version 3.1.0 (2014-04-10)
Platform: x86_64-apple-darwin13.1.0 (64-bit)
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] ggplot2_1.0.0 plyr_1.8.1 rlecuyer_0.3-3 snow_0.3-13
[5] qtl_1.32-10 knitr_1.6
loaded via a namespace (and not attached):
[1] colorspace_1.2-4 digest_0.6.4 evaluate_0.5.5 formatR_0.10
[5] grid_3.1.0 gtable_0.1.2 MASS_7.3-33 munsell_0.4.2
[9] parallel_3.1.0 proto_0.3-10 Rcpp_0.11.2 reshape2_1.4
[13] scales_0.2.4 stringr_0.6.2 tools_3.1.0
#############################
###Working Directory###
#############################
rm(list=ls())
setwd("~/Desktop/Drydown2014")
cross<-read.cross("csvs",dir="",
genotypes=c("a","b"),
genfile="TKrils_168_jm_gen.csv",
phefile= "TKrils_full_phe.csv",
na.strings=c("NA","-", "#VALUE!","#NUM!"))
--Read the following data:
341 individuals
168 markers
38 phenotypes
--Cross type: bc
#1.2: Set class as recombinant inbred line (DH is not an option in R/qtl)
class(cross)[1] <-"riself"
summary(cross)
RI strains via selfing
No. individuals: 341
No. phenotypes: 38
Percent phenotyped: 78.7
No. chromosomes: 5
Autosomes: 1 2 3 4 5
Total markers: 168
No. markers: 41 33 26 32 36
Percent genotyped: 91.1
Genotypes (%): AA:50.8 BB:49.2
#1.3 cull matrix to include only high quality markers
#look at the genotypes... check for errors
par(mfrow=c(1,5),mar = c(2,0,2,0))
cross.errlod<-calc.errorlod(cross)
for (i in 1:5) plotGeno(cross.errlod,chr=i, include.xo=F, cutoff=4, main=paste("ch",i))
# qtl analysis of regular map indicates that there is an error on ch4, between 45 and 46
print(toperr <- top.errorlod(cross.errlod, cutoff=4))
chr id marker errorlod
1 5 182 5_4223381 4.833
2 5 194 5_4223381 4.669
3 5 154 MSAT5.14 4.387
4 5 133 MSAT5.14 4.387
5 4 108 MSAT4.35 4.211
6 2 55 2_4354761 4.165
cross.clean <- cross
for(i in 1:nrow(toperr)) {
chr <- toperr$chr[i]; id <- toperr$id[i]; mar <- toperr$marker[i]
cross.clean$geno[[chr]]$data[cross$pheno$id==id, mar] <- NA
}
for (i in 1:5) plotGeno(calc.errorlod(cross.clean),chr=i, include.xo=T,cutoff=4, main=paste("ch",i))
#if this looks good, replace the old cross
##############
#still problems with chr 4 44-48#
cross<-cross.clean
#make pseudomarker matrix
grid<-calc.genoprob(cross, step=1, error.prob=.0001, map.function="kosambi",
stepwidth="max", off.end=0)
mainphes<-names(cross$phe)[-c(1,seq(from=4,to=28,by=3),33,36)] #all environment-specific QTLs
#make cross/pseudomarker matrix with quantile normalized phenotypes
quantnorm<-function(x) {
n=sum(!is.na(x),na.rm=T)
x=rank(x)/(n+1)
x=qnorm(x)
x[is.infinite(x)]=NA
x
}
cross.qn<-transformPheno(cross,pheno.col=mainphes,transf=quantnorm)
grid.qn<-calc.genoprob(cross.qn, step=1, error.prob=.0001, map.function="kosambi",
stepwidth="max", off.end=0)
#compare permutation results between groups
permtest.raw<-scanone(grid, method="hk",pheno.col=mainphes,n.perm=1000)
Permutation 20
Permutation 40
Permutation 60
Permutation 80
Permutation 100
Permutation 120
Permutation 140
Permutation 160
Permutation 180
Permutation 200
Permutation 220
Permutation 240
Permutation 260
Permutation 280
Permutation 300
Permutation 320
Permutation 340
Permutation 360
Permutation 380
Permutation 400
Permutation 420
Permutation 440
Permutation 460
Permutation 480
Permutation 500
Permutation 520
Permutation 540
Permutation 560
Permutation 580
Permutation 600
Permutation 620
Permutation 640
Permutation 660
Permutation 680
Permutation 700
Permutation 720
Permutation 740
Permutation 760
Permutation 780
Permutation 800
Permutation 820
Permutation 840
Permutation 860
Permutation 880
Permutation 900
Permutation 920
Permutation 940
Permutation 960
Permutation 980
Permutation 1000
permtest.qn<-scanone(grid.qn, method="hk",pheno.col=mainphes,n.perm=1000)
Permutation 20
Permutation 40
Permutation 60
Permutation 80
Permutation 100
Permutation 120
Permutation 140
Permutation 160
Permutation 180
Permutation 200
Permutation 220
Permutation 240
Permutation 260
Permutation 280
Permutation 300
Permutation 320
Permutation 340
Permutation 360
Permutation 380
Permutation 400
Permutation 420
Permutation 440
Permutation 460
Permutation 480
Permutation 500
Permutation 520
Permutation 540
Permutation 560
Permutation 580
Permutation 600
Permutation 620
Permutation 640
Permutation 660
Permutation 680
Permutation 700
Permutation 720
Permutation 740
Permutation 760
Permutation 780
Permutation 800
Permutation 820
Permutation 840
Permutation 860
Permutation 880
Permutation 900
Permutation 920
Permutation 940
Permutation 960
Permutation 980
Permutation 1000
permsout.raw<-data.frame(phe=character(),perm.raw.1=numeric(),perm.raw.05=numeric())
for (i in 1:length(mainphes)){
phe<-mainphes[i]
index.05<-(1:length(mainphes)*2)[i]-1
index.1<-(1:length(mainphes)*2)[i]
perm.raw.1<-summary(permtest.raw)[index.1]
perm.raw.05<-summary(permtest.raw)[index.05]
permsout.raw<-rbind(permsout.raw,data.frame(phe,perm.raw.1,perm.raw.05))
}
permsout.qn<-data.frame(phe=character(),perm.qn.1=numeric(),perm.qn.05=numeric())
for (i in 1:length(mainphes)){
phe<-mainphes[i]
index.05<-(1:length(mainphes)*2)[i]-1
index.1<-(1:length(mainphes)*2)[i]
perm.qn.1<-summary(permtest.qn)[index.1]
perm.qn.05<-summary(permtest.qn)[index.05]
permsout.qn<-rbind(permsout.qn,data.frame(phe,perm.qn.1,perm.qn.05))
}
outmaxlod<-data.frame(phe=character(),maxlod_raw=numeric(),maxlod_qn=numeric())
for (i in mainphes){
phe<-i
maxlod.raw<-max(scanone(grid, pheno.col=i, method="hk"))[1,3]
maxlod.qn<-max(scanone(grid.qn, pheno.col=i, method="hk"))[1,3]
outmaxlod<-rbind(outmaxlod,data.frame(phe,maxlod.raw,maxlod.qn))
}
anyqtls<-merge(permsout.raw,permsout.qn,by="phe")
anyqtls<-merge(anyqtls,outmaxlod,by="phe")
print(anyqtls) # the distribution of max scanone peaks and permutations
phe perm.raw.1 perm.raw.05 perm.qn.1 perm.qn.05 maxlod.raw
1 ABAlyoph.dry 2.227 2.490 2.325 2.641 3.5237
2 ABAlyoph.wet 2.210 2.512 2.261 2.539 1.9033
3 C13.dry 2.384 2.765 2.324 2.698 2.9466
4 C13.wet 2.363 2.635 2.229 2.427 5.5038
5 FT.June 2.318 2.628 2.271 2.624 53.0746
6 g0 2.341 2.651 2.314 2.651 4.4961
7 GR.dry 2.441 2.672 2.277 2.567 3.7482
8 LeafArea.dry 2.235 2.589 2.320 2.607 2.4588
9 LeafArea.wet 2.281 2.570 2.252 2.544 2.7781
10 proline.dry 2.281 2.576 2.307 2.621 16.0433
11 proline.wet 2.215 2.507 2.318 2.638 1.5173
12 Rel.GR.wet 2.314 2.639 2.277 2.603 3.4307
13 RMR.dry 2.427 2.794 2.488 2.800 2.6379
14 RMR.wet 2.413 2.710 2.376 2.670 1.5504
15 RtMass.dry 2.438 2.786 2.360 2.767 2.3300
16 RtMass.wet 2.311 2.640 2.410 2.782 1.3313
17 S.R.dry 2.468 2.764 2.440 2.810 1.5904
18 S.R.wet 2.162 2.479 2.405 2.745 1.3005
19 ShtDW.dry 2.366 2.768 2.319 2.577 1.7355
20 ShtDW.wet 2.244 2.528 2.290 2.659 1.5474
21 ShtFrWt.dry 2.286 2.619 2.301 2.568 2.1435
22 ShtFrWt.wet 2.327 2.634 2.320 2.600 1.8202
23 WC.dry 2.220 2.525 2.366 2.606 2.3468
24 WC.wet 2.317 2.571 2.342 2.711 0.9362
25 Wilted.dry 2.276 2.667 2.300 2.530 2.7511
26 Wilted.wet 2.036 2.266 2.023 2.232 1.2578
maxlod.qn
1 2.4009
2 1.4075
3 2.9618
4 5.6963
5 53.7338
6 4.5293
7 3.5959
8 2.4290
9 2.6025
10 16.4985
11 2.2366
12 2.8481
13 2.3096
14 1.3736
15 2.2925
16 0.9848
17 1.5294
18 1.1509
19 1.6709
20 1.4349
21 2.2253
22 1.9064
23 2.2975
24 1.0334
25 2.7244
26 1.2578
anyqtls$qtl.qn.05<-anyqtls$maxlod.qn>=anyqtls$perm.qn.05
anyqtls$qtl.raw.05<-anyqtls$maxlod.raw>=anyqtls$perm.raw.05
anyqtls$qtl.qn.1<-anyqtls$maxlod.qn>=anyqtls$perm.qn.1
anyqtls$qtl.raw.1<-anyqtls$maxlod.raw>=anyqtls$perm.raw.1
print(anyqtls[c(1,8:11)])
phe qtl.qn.05 qtl.raw.05 qtl.qn.1 qtl.raw.1
1 ABAlyoph.dry FALSE TRUE TRUE TRUE
2 ABAlyoph.wet FALSE FALSE FALSE FALSE
3 C13.dry TRUE TRUE TRUE TRUE
4 C13.wet TRUE TRUE TRUE TRUE
5 FT.June TRUE TRUE TRUE TRUE
6 g0 TRUE TRUE TRUE TRUE
7 GR.dry TRUE TRUE TRUE TRUE
8 LeafArea.dry FALSE FALSE TRUE TRUE
9 LeafArea.wet TRUE TRUE TRUE TRUE
10 proline.dry TRUE TRUE TRUE TRUE
11 proline.wet FALSE FALSE FALSE FALSE
12 Rel.GR.wet TRUE TRUE TRUE TRUE
13 RMR.dry FALSE FALSE FALSE TRUE
14 RMR.wet FALSE FALSE FALSE FALSE
15 RtMass.dry FALSE FALSE FALSE FALSE
16 RtMass.wet FALSE FALSE FALSE FALSE
17 S.R.dry FALSE FALSE FALSE FALSE
18 S.R.wet FALSE FALSE FALSE FALSE
19 ShtDW.dry FALSE FALSE FALSE FALSE
20 ShtDW.wet FALSE FALSE FALSE FALSE
21 ShtFrWt.dry FALSE FALSE FALSE FALSE
22 ShtFrWt.wet FALSE FALSE FALSE FALSE
23 WC.dry FALSE FALSE FALSE TRUE
24 WC.wet FALSE FALSE FALSE FALSE
25 Wilted.dry TRUE TRUE TRUE TRUE
26 Wilted.wet FALSE FALSE FALSE FALSE
phes_wqtls<-as.character(anyqtls$phe[anyqtls$qtl.raw.05])
phes_wqtls #these are the traits with qtl at a=0.05 w/ raw phenotypes
[1] "ABAlyoph.dry" "C13.dry" "C13.wet" "FT.June"
[5] "g0" "GR.dry" "LeafArea.wet" "proline.dry"
[9] "Rel.GR.wet" "Wilted.dry"
#compare scanone of phenotypes between qn and raw
par(mfrow=c(3,1),mar = c(2,4,2,0))
for (i in phes_wqtls){
s1.raw<-scanone(grid, pheno.col=i)
s1.qn<-scanone(grid.qn, pheno.col=i)
plot(s1.raw, ylim=c(0,max(max(s1.qn$lod),(.5+max(s1.raw$lod)))),
main=paste(i, " red: quantile normalized"), mar=c(0,0,0,0))
plot(s1.qn, col="red", lty=2, add=T)
}
min.thresh<-min(anyqtls$perm.raw.05) #most lenient QTL selection criteria
# Part 2:
# heavy and light penalties from previous analyses
pens<-c(min.thresh,3.2305306,1.68988784)
#qtlphes<-as.character(outmaxlod[outmaxlod$maxlod>=pens[1],]$i)
qtlphes<-phes_wqtls #these are the traits w/qtls
#set up the environment to write statistics to
stats_out<-data.frame(phenotype=character(),chromosome=numeric(),position=numeric(),
df=numeric(),type3SS=numeric(),LOD=numeric(),perc.var=numeric(),Fstat=numeric(),P.chi2=numeric(),P.F=numeric(),
effect.estimate=numeric(),effect.SE=numeric(),effect.t=numeric(),
lowCImarker=character(), hiCImarker=character(),
lowCIpos=character(), hiCIpos=character())
statcols<-names(stats_out)
par(mfrow=c(4,1),mar = c(2,4,2,1))
#this is the engine that generates statistics for each phenotype
for (i in qtlphes){
print(i)
stepout<-stepwiseqtl(grid,pheno.col=i, penalties=pens, max.qtl=10,method="hk",
refine.locations=T,additive.only=T, keeplodprofile=T,keeptrace=T,verbose=F)
plotLodProfile(stepout,main=i)
fit<-fitqtl(grid,pheno.col=i,qtl=stepout,formula=formula(stepout),get.ests=T,dropone=T)
nqtls<-nqtl(stepout)
if(nqtls>1){
ests_out<-summary(fit)$ests[2:(nqtls+1),1:3]
cis<-data.frame()
for (j in 1:nqtls){
ciout<-bayesint(stepout,qtl.index=j, expandtomarkers=T, prob=.9)
lowmarker<-rownames(ciout)[1]; highmarker<-rownames(ciout)[3]
lowposition<-ciout[1,2]; highposition<-ciout[3,2]
cis<-rbind(cis,cbind(lowmarker,highmarker,lowposition,highposition))
}
stats<-data.frame(rep(i,nqtls),stepout$chr,stepout$pos,fit$result.drop,ests_out,cis)
colnames(stats)<-statcols
}
else{
cis<-data.frame()
for (j in 1:nqtls){
ciout<-bayesint(stepout,qtl.index=j, expandtomarkers=T, prob=.9)
lowmarker<-rownames(ciout)[1]; highmarker<-rownames(ciout)[3]
lowposition<-ciout[1,2]; highposition<-ciout[3,2]
cis<-rbind(cis,cbind(lowmarker,highmarker,lowposition,highposition))
}
ests_out<-summary(fit)$ests[2:(nqtls+1),1:3]
stats<-data.frame(c(i,stepout$chr[1],stepout$pos[1],as.numeric(fit$result.full[1,c(1,2,4,5,3,6,7)]),as.numeric(ests_out),cis[1,]))
colnames(stats)<-statcols; rownames(stats)<-stepout$name
}
stats_out<-rbind(stats_out,stats)
}
[1] "ABAlyoph.dry"
[1] "C13.dry"
[1] "C13.wet"
[1] "FT.June"
[1] "g0"
[1] "GR.dry"
[1] "LeafArea.wet"
[1] "proline.dry"
[1] "Rel.GR.wet"
[1] "Wilted.dry"
write.csv(stats_out, file="dd2014_allqtlstatistics.csv")
print(stats_out) # this is the statistical output in a dataframe
phenotype chromosome position df type3SS LOD perc.var
2@13.5 ABAlyoph.dry 2 13.469 1 2.620e+02 3.524 5.006
2@70.3 C13.dry 2 70.323 1 4.412e+00 3.463 4.265
3@47.5 C13.dry 3 47.460 1 4.554e+00 3.572 4.402
4@37.0 C13.dry 4 37.031 1 3.426e+00 2.703 3.311
5@72.2 C13.dry 5 72.223 1 4.732e+00 3.708 4.574
1@60.8 C13.wet 1 60.778 1 2.074e+00 3.009 3.029
3@15.6 C13.wet 3 15.551 1 1.726e+00 2.514 2.522
3@90.2 C13.wet 3 90.218 1 2.630e+00 3.795 3.842
4@1.0 C13.wet 4 0.951 1 5.761e+00 8.069 8.417
4@45.1 C13.wet 4 45.123 1 2.723e+00 3.927 3.979
5@11.8 C13.wet 5 11.761 1 3.396e+00 4.865 4.962
5@73.2 C13.wet 5 73.162 1 4.576e+00 6.481 6.685
1@30.6 FT.June 1 30.579 1 4.181e+02 13.534 5.905
1@31.5 FT.June 1 31.547 1 3.053e+02 10.121 4.312
1@83.6 FT.June 1 83.628 1 3.535e+02 11.599 4.993
4@0.0 FT.June 4 0.000 1 4.350e+03 83.286 61.440
4@65.9 FT.June 4 65.929 1 2.112e+02 7.149 2.984
5@15.2 FT.June 5 15.168 1 1.075e+02 3.726 1.519
5@69.4 FT.June 5 69.411 1 3.423e+02 11.258 4.835
1@0.0 g0 1 0.000 1 6.189e-03 5.476 7.789
1@99.3 g0 1 99.350 1 2.968e-03 2.686 3.735
2@31.3 g0 2 31.312 1 4.062e-03 3.648 5.112
3@25.6 GR.dry 3 25.589 1 1.436e+02 3.748 5.283
4@45.11 LeafArea.wet 4 45.123 1 6.686e+02 2.528 3.264
4@46.0 LeafArea.wet 4 46.041 1 1.146e+03 4.282 5.594
2@30.4 proline.dry 2 30.360 1 3.268e+04 2.518 2.472
2@65.6 proline.dry 2 65.575 1 2.765e+05 18.903 20.922
3@63.4 proline.dry 3 63.395 1 5.096e+04 3.888 3.855
4@1.9 proline.dry 4 1.902 1 4.393e+04 3.365 3.324
4@81.2 proline.dry 4 81.217 1 3.175e+04 2.448 2.402
3@86.6 Rel.GR.wet 3 86.605 1 3.776e-02 3.431 4.862
2@14.3 Wilted.dry 2 14.325 1 6.583e+01 2.885 3.721
4@43.4 Wilted.dry 4 43.385 1 9.915e+01 4.301 5.605
4@48.0 Wilted.dry 4 47.996 1 1.411e+02 6.045 7.979
Fstat P.chi2 P.F effect.estimate effect.SE effect.t
2@13.5 261.96100 5.618e-05 6.011e-05 -0.982492 0.241537 -4.068
2@70.3 16.09396 6.510e-05 7.505e-05 0.119176 0.029707 4.012
3@47.5 16.61239 4.997e-05 5.785e-05 -0.120547 0.029576 -4.076
4@37.0 12.49585 4.181e-04 4.678e-04 -0.104069 0.029440 -3.535
5@72.2 17.25985 3.596e-05 4.185e-05 0.137703 0.033145 4.154
1@60.8 13.80889 1.973e-04 2.381e-04 0.088624 0.023849 3.716
3@15.6 11.49692 6.680e-04 7.828e-04 0.076349 0.022517 3.391
3@90.2 17.51141 2.910e-05 3.677e-05 -0.090149 0.021543 -4.185
4@1.0 38.36420 1.088e-09 1.767e-09 -0.137756 0.022241 -6.194
4@45.1 18.13610 2.115e-05 2.694e-05 -0.094149 0.022108 -4.259
5@11.8 22.61612 2.211e-06 2.975e-06 -0.106421 0.022378 -4.756
5@73.2 30.47367 4.675e-08 6.916e-08 0.136853 0.024791 5.520
1@30.6 66.80920 2.887e-15 6.439e-15 3.539286 0.433010 8.174
1@31.5 48.78444 8.673e-12 1.574e-11 -2.935921 0.420343 -6.985
1@83.6 56.48653 2.702e-13 5.340e-13 1.089161 0.144917 7.516
4@0.0 695.10692 0.000e+00 0.000e+00 -3.669350 0.139176 -26.365
4@65.9 33.75760 9.594e-09 1.467e-08 -0.846001 0.145608 -5.810
5@15.2 17.18333 3.441e-05 4.316e-05 -0.582742 0.140580 -4.145
5@69.4 54.69717 6.007e-13 1.164e-12 1.098149 0.148484 7.396
1@0.0 25.99308 5.117e-07 6.280e-07 0.004986 0.000978 5.098
1@99.3 12.46494 4.361e-04 4.839e-04 0.003492 0.000989 3.531
2@31.3 17.06047 4.153e-05 4.771e-05 -0.004837 0.001171 -4.130
3@25.6 143.60123 3.258e-05 3.499e-05 -0.754263 0.179655 -4.198
4@45.11 11.74181 6.441e-04 6.872e-04 -4.037755 1.178345 -3.427
4@46.0 20.12427 8.975e-06 9.987e-06 5.242976 1.168740 4.486
2@30.4 11.58674 6.612e-04 7.496e-04 -13.499903 3.965976 -3.404
2@65.6 98.05592 0.000e+00 0.000e+00 33.785459 3.411873 9.902
3@63.4 18.06865 2.322e-05 2.806e-05 -14.560774 3.425481 -4.251
4@1.9 15.57780 8.273e-05 9.756e-05 -11.797479 2.989070 -3.947
4@81.2 11.25811 7.868e-04 8.891e-04 -10.008236 2.982804 -3.355
3@86.6 0.03776 7.043e-05 7.522e-05 -0.011273 0.002810 -4.012
2@14.3 13.39481 2.675e-04 2.955e-04 0.502782 0.137376 3.660
4@43.4 20.17576 8.569e-06 9.908e-06 -1.281149 0.285223 -4.492
4@48.0 28.71884 1.319e-07 1.613e-07 1.540130 0.287392 5.359
lowCImarker hiCImarker lowCIpos hiCIpos
2@13.5 2_518639 2_9256494 7.378 40.829
2@70.3 2-15986145 2_17844202 67.111 72.187
3@47.5 MSAT3.5 3.F1P2-TGF 40.691 68.384
4@37.0 4_4074232 4-17608717 26.868 81.217
5@72.2 5-22414510 5_26121185 68.476 83.371
1@60.8 1-17355484 1.NGA128 52.749 70.616
3@15.6 3.nga162 MSAT3.6 0 38.289
3@90.2 3_21728447 3-20482976 87.311 94.929
4@1.0 4.FRI 4_2126308 0 24.714
4@45.1 MSAT4.16 4-12403706 40.399 57.095
5@11.8 5.MUK11D 5.NGA106 6.776 15.168
5@73.2 5-22414510 5_25301310 68.476 79.489
1@30.6 1-9343234 1-9343234 30.579 30.579
1@31.5 1-9496666 1-9496666 31.547 31.547
1@83.6 1.NGA128 MSAT1.1 70.616 92.834
4@0.0 4.FRI 4.FRI 0 0
4@65.9 4_12900428 MSAT4.30 59.969 74.562
5@15.2 5_2799359 5-7278798 10.909 25.597
5@69.4 5-22414510 5-25301036 68.476 75.978
1@0.0 1_444764 1_112908 0 3.027
1@99.3 1.F5I14 1.NGA692 82.707 104.043
2@31.3 2_518639 2_9256494 7.378 40.829
3@25.6 3_620734 3-7391896 6.849 31.898
4@45.11 4_5726702 4-8078388 31.979 45.123
4@46.0 4-8297594 4-10482038 46.041 49.952
2@30.4 2-2477306 2_9256494 18.081 40.829
2@65.6 2_14801461 2-15986145 64.807 67.111
3@63.4 3.GAPABX 3.F1P2-TGF 47.46 68.384
4@1.9 4.FRI 4_2126308 0 24.714
4@81.2 4_208623 4-17608717 1.902 81.217
3@86.6 3-9747158 3-21632379 51.096 96.869
2@14.3 2_518639 2_9256494 7.378 40.829
4@43.4 MSAT4.16 4-8078388 40.399 45.123
4@48.0 4-8297594 4-10482038 46.041 49.952
#Part 2: Processing and plotting the QTL statistics
#2.1 Read in all relevant data
statsout<-read.csv("dd2014_allqtlstatistics.csv",header=T) #statistics from Part 1
#2.2 Plot some phenotypes, determine correlations between phenotypes
phematrix<-cross$phe[,qtlphes]
distmatrix<-dist(t(phematrix))
hc<-hclust(distmatrix)
par(mfrow=c(1,1))
plot(hc)#plot the similarity of phenotypes
#reorder the phenotypes so that they are clustered by distance
statsout$phenotype <- factor(statsout$phenotype,
levels = names(phematrix)[hc$order])
#2.3plot the distribution of QTL
#
ggplot(statsout,aes(x=position,y=as.numeric(effect.t)))+
geom_point() +
scale_y_continuous("t statistic for QTL point estimate")+
geom_errorbarh(aes(xmax = hiCIpos, xmin = lowCIpos, height = 0))+
facet_grid(phenotype~chromosome, scales="free_y", margins=F)+
theme_bw()+
ggtitle("distribution of QTL positions and effect sizes")
ggplot(statsout,aes(x=position,y=P.F))+
geom_point() +
scale_y_log10("p.value of F statistic")+
geom_errorbarh(aes(xmax = hiCIpos, xmin = lowCIpos, height = 0))+
facet_grid(.~chromosome, scales="free_y", margins=F)+
theme_bw()+
ggtitle("distribution of QTL positions and effect sizes")
ggplot(statsout,aes(x=position))+
geom_histogram(binwidth=20) +
facet_grid(.~chromosome, scales="free_y", margins=F)+
theme_bw()+
ggtitle("distribution and density of QTL point estimates")
