Lovell — Aug 11, 2014, 6:49 PM
#determination of QTL intervals for candidate gene discovery
rm(list=ls())
options(warn=-1)
pkg <- c("biomaRt","qtl","snow","rlecuyer","plyr","ggplot2","reshape","reshape2","qvalue","lme4","nlme","car","doBy","scales")
invisible(lapply(pkg, function(x) {cat(x,"..."); library(x, character.only=T, verbose=F, warn.conflicts=F,quietly=T)} ))
biomaRt ...qtl ...snow ...rlecuyer ...plyr ...ggplot2 ...reshape ...reshape2 ...qvalue ...lme4 ...
Attaching package: 'Matrix'
The following object is masked from 'package:reshape':
expand
nlme ...car ...doBy ...
Loading required package: splines
Attaching package: 'MASS'
The following object is masked from 'package:biomaRt':
select
scales ...
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] splines stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] scales_0.2.4 doBy_4.5-10 MASS_7.3-33 survival_2.37-7
[5] car_2.0-20 nlme_3.1-117 lme4_1.1-7 Rcpp_0.11.2
[9] Matrix_1.1-4 qvalue_1.38.0 reshape2_1.4 reshape_0.8.5
[13] ggplot2_1.0.0 plyr_1.8.1 rlecuyer_0.3-3 snow_0.3-13
[17] qtl_1.32-10 biomaRt_2.20.0 knitr_1.6
loaded via a namespace (and not attached):
[1] AnnotationDbi_1.26.0 Biobase_2.24.0 BiocGenerics_0.10.0
[4] colorspace_1.2-4 DBI_0.2-7 digest_0.6.4
[7] evaluate_0.5.5 formatR_0.10 GenomeInfoDb_1.0.2
[10] grid_3.1.0 gtable_0.1.2 IRanges_1.22.10
[13] lattice_0.20-29 minqa_1.2.3 munsell_0.4.2
[16] nloptr_1.0.0 nnet_7.3-8 parallel_3.1.0
[19] proto_0.3-10 RCurl_1.95-4.3 RSQLite_0.11.4
[22] stats4_3.1.0 stringr_0.6.2 tcltk_3.1.0
[25] tools_3.1.0 XML_3.98-1.1
setwd("~/Desktop/Drydown2014")
##############################################
##############################################
#Part 1 (QTL analysis w/out covariates):
##############################################
##############################################
##############################################
load("dd2014_stepwiseresults.RData")
library(qtl)
library(ggplot2)
library (plyr)
#1.1: Look at the distribution of QTL confidence intervals
plod.nocovar<-lapply(model.out, function(x) attr(x, "pLOD"))
plod.cytocovar<-lapply(model.out.cyto, function(x) attr(x, "pLOD"))
plods<-cbind(unlist(plod.cytocovar),unlist(plod.nocovar)); colnames(plods)<-c("cytocovar","nocovar")
cyto.is.best<-plods[,1]>plods[,2]; cyto.is.best["FT.June"]<-F
cyto.is.best<-as.character(names(cyto.is.best[cyto.is.best]))
cyto.models<-model.out.cyto[cyto.is.best]
nocov.is.best<-plods[,1]<=plods[,2];nocov.is.best["FT.June"]<-T
nocov.is.best<-as.character(names(nocov.is.best[nocov.is.best]))
no.models<-model.out[nocov.is.best]
bestmodels<-append(cyto.models,no.models)
#extract all lps into a dataframe
all.lps<-lapply(bestmodels, function (x) attr (x, "lodprofile"))
names(all.lps)
[1] "GR.dry" "Rel.GR.wet" "Rolled.dry" "Wilted.dry"
[5] "g0" "proline.dry" "ABAlyoph.dry" "C13.dry"
[9] "C13.wet" "FT.June" "RMR.dry"
lps.df<-data.frame()
count.1<-0
for (i in all.lps){
count.1=count.1+1
count=0
for (j in i){
count=count+1
out<-cbind(names(all.lps)[count.1],names(i)[count],as.data.frame(j))
st.lod<-out$lod/max(out$lod)
out.all<-cbind(out,st.lod)
lps.df<-rbind(lps.df,out.all)
}
}
colnames(lps.df)<-c("phenotype","qtl.id","chr","pos","lod","st.lod")
lps.df$lod[lps.df$lod<0]<-NA
lps.df$st.lod[lps.df$st.lod<0]<-NA
# extract statistics for each model and combine with Lod Profiles
cyto.is.best;nocov.is.best
[1] "GR.dry" "Rel.GR.wet" "Rolled.dry" "Wilted.dry" "g0"
[6] "proline.dry"
[1] "ABAlyoph.dry" "C13.dry" "C13.wet" "FT.June"
[5] "RMR.dry"
plod.cytocovar<-lapply(cytocovar.modelstats, function(x) attr(x, "pLOD"))
plod.nocovar<-lapply(nocovar.modelstats, function(x) attr(x, "pLOD"))
cyto.stats.in<-split(cytocovar.modelstats,cytocovar.modelstats$phenotype)
nocov.stats.in<-split(nocovar.modelstats,nocovar.modelstats$phenotype)
cyto.stats<-cyto.stats.in[cyto.is.best]
nocov.stats<-nocov.stats.in[nocov.is.best]
beststats<-append(cyto.stats,nocov.stats)
beststats.df<-ldply(beststats,data.frame)
beststats.ci<-beststats.df[complete.cases(beststats.df),c("phenotype","chromosome","position",
"lowCImarker", "hiCImarker", "lowCIpos", "hiCIpos")]
#get qtl identifiers and add them to the data frame
qtlids<-data.frame()
for (i in 1:11){
temp<-bestmodels[[i]]; qtl.id<-temp$name; qtl.pos<-temp$pos; qtl.chr<-temp$chr
pheno<-names(bestmodels)[i]
qtlids<-rbind(qtlids,cbind(qtl.id,qtl.chr,qtl.pos,pheno))
}
names(qtlids)[2:4]<-c("chromosome","position","phenotype")
allstats<-merge(qtlids,beststats.ci, by=c("chromosome","position","phenotype"))
all.stats.out<-merge(lps.df,allstats, by=c("qtl.id","phenotype"), all.x=T)
all.stats.out$unique<-paste(all.stats.out$qtl.id,all.stats.out$phenotype, sep="_")
test.oneelement<-all.stats.out[!duplicated(all.stats.out$unique),]
write.csv(test.oneelement, file="dd2014_ciintervals.csv", row.names=F)
write.csv(all.stats.out, file="dd2014_ciandlps.csv", row.names=F)
for.rects<-read.csv("dd2014_ciintervals.csv", header=T)
for.lps<-read.csv("dd2014_ciandlps.csv", header=T)
ggplot(for.rects, aes(x=pos, y=st.lod, fill=qtl.id, color=qtl.id))+
geom_rect(aes(xmax=hiCIpos, xmin=lowCIpos, ymin = 0, ymax = 1.02), color=NA,fill="red", alpha=.2)+
geom_line(aes(x=pos, y=st.lod),color="black",data=for.lps)+
facet_grid(phenotype~chromosome, space = "free", scales= "free")+
theme_bw()+
theme(
plot.background = element_blank()
,panel.grid.major = element_blank()
,panel.grid.minor = element_blank(),
strip.background = element_rect(colour="white", fill="white"),
panel.border = element_rect(colour = "black"))+
scale_y_continuous(breaks=c(0,.5,1))+
ggtitle("distribution of standardized LOD profiles overlayed by 90% CIs")
## now look at specific problems
print(for.rects$unique)
[1] 1@16.1_g0 1@83.8_FT.June 2@15.5_ABAlyoph.dry
[4] 2@73.6_Rolled.dry 2@74.4_proline.dry 3@17.9_C13.wet
[7] 3@19.0_RMR.dry 3@3.9_GR.dry 3@50.2_proline.dry
[10] 3@53.6_Rel.GR.wet 3@58.5_C13.dry 3@75.5_C13.wet
[13] 4@1.5_RMR.dry 4@2.1_Wilted.dry 4@2.8_FT.June
[16] 4@2.8_proline.dry 4@4.4_C13.wet 4@42.4_C13.wet
[19] 4@62.1_FT.June 5@15.3_FT.June 5@69.8_C13.wet
[22] 5@71.7_FT.June 5@9.8_C13.wet
23 Levels: 1@16.1_g0 1@83.8_FT.June ... 5@9.8_C13.wet
#make lists- qtl.index, phenotype
too.wide<-list(c(1,"ABAlyoph.dry"),c(1,"GR.dry"))
too.narrow<-list(c(1,"proline.dry"),c(2,"FT.June"),c(5,"FT.June"))
too.wide.unique<-c("2@15.5_ABAlyoph.dry","3@3.9_GR.dry")
too.narrow.unique<-c("2@74.4_proline.dry","4@2.8_FT.June","5@71.7_FT.June")
newint<-data.frame(bayesint(model.out["ABAlyoph.dry"][[1]], qtl.index=1, prob=.75)[1:3,1:3])
#function to look at plots of new vs old cis.
fix.qtlci<-function(qtllist,qtlnames, newdrop=1.5,newprob=.90) {
all.out<-data.frame()
for (i in 1:length(qtllist)){
for.rects.sub<-for.rects[for.rects$unique==qtlnames[i],]
for.lps.sub<-for.lps[for.lps$unique==qtlnames[i],]
newbayesint<-data.frame(bayesint(bestmodels[qtllist[[i]][2]][[1]],
qtl.index=as.numeric(qtllist[[i]][1]), prob=newprob))
newlodint<-data.frame(lodint(bestmodels[qtllist[[i]][2]][[1]],
qtl.index=as.numeric(qtllist[[i]][1]), drop=newdrop))
print(ggplot(for.rects.sub, aes(x=pos, y=st.lod, fill=qtl.id, color=qtl.id))+
geom_rect(aes(xmax=hiCIpos, xmin=lowCIpos, ymin = 0, ymax = 1.02), color="green",fill=NA)+
geom_line(aes(x=pos, y=st.lod),color="black",data=for.lps.sub)+theme_bw()+
geom_rect(xmax=newlodint[3,2],xmin=newlodint[1,2],ymin=0,ymax=1.02,
fill="blue",color=NA, alpha=0.2)+
geom_rect(xmax=newbayesint[3,2],xmin=newbayesint[1,2],ymin=0,ymax=1.02,
fill="red",color=NA, alpha=0.2)+
ggtitle(paste(qtlnames[i], "original 90% CI (green box) vs new CIs \n", newprob*100 ,"% bayes (red) and", newdrop, "drop (blue)"))
)
ints.out<-cbind(strsplit(qtlnames[i],"_")[[1]][1],
strsplit(qtlnames[i],"_")[[1]][2],
newlodint,newbayesint)
all.out<-rbind(all.out,ints.out)
}
colnames(all.out)[1:2]<-c("qtl.id","phenotype")
colnames(all.out)[3:5]<-paste(colnames(all.out)[3:5],"_lod",sep="")
colnames(all.out)[6:8]<-paste(colnames(all.out)[6:8],"_bayes",sep="")
return(all.out)
}
#this procedure makes estimates of extremely narrow peaks more reasonable and
# and broad peaks more representative of the morphology
toowide.newcis<-fix.qtlci(qtllist=too.wide, qtlnames=too.wide.unique, newdrop=1,newprob=.75)
toonarrow.newcis<-fix.qtlci(qtllist=too.narrow, qtlnames=too.narrow.unique, newdrop=2,newprob=.99)
#best ints
inttoo.wide.unique<-c("2@15.5_ABAlyoph.dry","3@3.9_GR.dry")
too.narrow.unique<-c("2@74.4_proline.dry","4@2.8_FT.June","5@71.7_FT.June")
all.newcis<-rbind(toowide.newcis,toonarrow.newcis)
Ft5.2<-"lod"
Ft4.1<-"lod"
pro2.2<-"bayes"
gr3.1<-"lod"
ABA2.1<-"bayes"
lod.newcis<-all.newcis[c(4:6,10:15),1:5]
bayes.newcis<-all.newcis[c(1:3,7:9),c(1:2,6:8)]
colnames(lod.newcis)[3:5]<-c("chr","pos","lod")
colnames(bayes.newcis)[3:5]<-c("chr","pos","lod")
##################
##################
#fixed cis:
newcis.out<-rbind(bayes.newcis,lod.newcis)
newcis.out$type<-rep(c("lowCIpos","position","hiCIpos"),5)
newcis.out<-dcast(newcis.out, qtl.id + phenotype +chr ~ type, value.var="pos")
##################
##################
#################
#now address the multi-peak problem
multi.peaks<-list(c(3,"C13.dry"),c(5,"C13.wet"),c(5,"C13.wet"),c(4,"Wilted.dry"), c(1,"g0"))
multi.peaks.unique<-c("3@58.5_C13.dry","5@9.8_C13.wet","5@69.8_C13.wet","4@2.1_Wilted.dry","1@16.1_g0")
#look at the standardized lps:
for (i in 1:5){
for.rects.sub<-for.rects[for.rects$unique==multi.peaks.unique[i],]
for.lps.sub<-for.lps[for.lps$unique==multi.peaks.unique[i],]
print(ggplot(for.rects.sub, aes(x=pos, y=st.lod, fill=qtl.id, color=qtl.id))+
geom_line(aes(x=pos, y=st.lod),color="black",data=for.lps.sub)+theme_bw()+
ggtitle(paste("lod profile for", multi.peaks.unique[i])))
}
#must do it manually- subset for peaks around 1/2 the height of the peak
new.ints=list()
#"3@58.5_C13.dry"
for.lps.sub<-for.lps[for.lps$unique=="3@58.5_C13.dry",]
new.int<-for.lps.sub[for.lps.sub$st.lod>=0.6,]
new.int
qtl.id phenotype chr pos lod st.lod chromosome position
1438 3@58.5 C13.dry 3 44.18 1.998 0.6178 3 58.47
1439 3@58.5 C13.dry 3 44.91 2.067 0.6391 3 58.47
1440 3@58.5 C13.dry 3 45.63 1.958 0.6054 3 58.47
1445 3@58.5 C13.dry 3 49.28 2.049 0.6335 3 58.47
1446 3@58.5 C13.dry 3 50.17 2.041 0.6312 3 58.47
1447 3@58.5 C13.dry 3 50.95 2.074 0.6414 3 58.47
1448 3@58.5 C13.dry 3 51.72 2.085 0.6446 3 58.47
1449 3@58.5 C13.dry 3 52.50 2.072 0.6408 3 58.47
1450 3@58.5 C13.dry 3 53.41 2.282 0.7057 3 58.47
1451 3@58.5 C13.dry 3 53.59 2.699 0.8346 3 58.47
1452 3@58.5 C13.dry 3 54.39 2.596 0.8027 3 58.47
1453 3@58.5 C13.dry 3 54.71 2.756 0.8522 3 58.47
1454 3@58.5 C13.dry 3 55.06 2.585 0.7993 3 58.47
1455 3@58.5 C13.dry 3 55.20 2.525 0.7806 3 58.47
1456 3@58.5 C13.dry 3 55.67 2.538 0.7849 3 58.47
1457 3@58.5 C13.dry 3 56.41 2.558 0.7910 3 58.47
1458 3@58.5 C13.dry 3 57.16 2.534 0.7836 3 58.47
1459 3@58.5 C13.dry 3 57.73 2.921 0.9033 3 58.47
1460 3@58.5 C13.dry 3 58.30 3.175 0.9818 3 58.47
1461 3@58.5 C13.dry 3 58.47 3.234 1.0000 3 58.47
1462 3@58.5 C13.dry 3 59.13 3.189 0.9859 3 58.47
1463 3@58.5 C13.dry 3 59.74 2.903 0.8978 3 58.47
1464 3@58.5 C13.dry 3 60.10 2.940 0.9092 3 58.47
1465 3@58.5 C13.dry 3 60.98 2.806 0.8676 3 58.47
1466 3@58.5 C13.dry 3 61.55 2.481 0.7672 3 58.47
1467 3@58.5 C13.dry 3 62.20 2.671 0.8258 3 58.47
1468 3@58.5 C13.dry 3 62.87 2.596 0.8026 3 58.47
1469 3@58.5 C13.dry 3 63.53 2.402 0.7426 3 58.47
1470 3@58.5 C13.dry 3 64.23 2.450 0.7575 3 58.47
1471 3@58.5 C13.dry 3 64.78 2.674 0.8270 3 58.47
1472 3@58.5 C13.dry 3 65.75 2.584 0.7991 3 58.47
1473 3@58.5 C13.dry 3 66.71 2.422 0.7489 3 58.47
1474 3@58.5 C13.dry 3 67.67 2.192 0.6779 3 58.47
1499 3@58.5 C13.dry 3 84.53 1.982 0.6129 3 58.47
1500 3@58.5 C13.dry 3 85.26 2.239 0.6923 3 58.47
1501 3@58.5 C13.dry 3 85.98 2.407 0.7444 3 58.47
1502 3@58.5 C13.dry 3 86.74 2.414 0.7464 3 58.47
1503 3@58.5 C13.dry 3 87.66 2.600 0.8039 3 58.47
1504 3@58.5 C13.dry 3 88.58 2.680 0.8286 3 58.47
1505 3@58.5 C13.dry 3 89.50 2.641 0.8166 3 58.47
lowCImarker hiCImarker lowCIpos hiCIpos unique
1438 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1439 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1440 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1445 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1446 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1447 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1448 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1449 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1450 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1451 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1452 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1453 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1454 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1455 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1456 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1457 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1458 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1459 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1460 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1461 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1462 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1463 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1464 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1465 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1466 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1467 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1468 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1469 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1470 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1471 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1472 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1473 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1474 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1499 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1500 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1501 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1502 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1503 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1504 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
1505 3.GAPABX 3-21632379 44.18 89.5 3@58.5_C13.dry
new.int<-new.int[new.int$pos> 49 & new.int$pos<68,]
new.ints[["3@58.5_C13.dry"]]<-c("3@58.5_C13.dry",min(new.int$pos),max(new.int$pos))
new.ints
$`3@58.5_C13.dry`
[1] "3@58.5_C13.dry" "49.2826666666667" "67.6726"
#"5@9.8_C13.wet"
for.lps.sub<-for.lps[for.lps$unique=="5@9.8_C13.wet",]
new.int<-for.lps.sub[for.lps.sub$st.lod>=0.6,]
new.int
qtl.id phenotype chr pos lod st.lod chromosome position
2783 5@9.8 C13.wet 5 5.785 3.065 0.6046 5 9.845
2784 5@9.8 C13.wet 5 6.442 3.145 0.6204 5 9.845
2786 5@9.8 C13.wet 5 7.410 3.195 0.6302 5 9.845
2787 5@9.8 C13.wet 5 8.095 3.563 0.7028 5 9.845
2788 5@9.8 C13.wet 5 8.779 3.716 0.7330 5 9.845
2789 5@9.8 C13.wet 5 9.606 4.617 0.9108 5 9.845
2790 5@9.8 C13.wet 5 9.845 5.070 1.0000 5 9.845
2791 5@9.8 C13.wet 5 10.396 5.002 0.9867 5 9.845
2792 5@9.8 C13.wet 5 10.948 4.527 0.8930 5 9.845
2793 5@9.8 C13.wet 5 11.575 4.373 0.8626 5 9.845
2794 5@9.8 C13.wet 5 11.641 4.367 0.8613 5 9.845
2795 5@9.8 C13.wet 5 11.955 4.552 0.8978 5 9.845
2796 5@9.8 C13.wet 5 12.797 4.053 0.7995 5 9.845
2797 5@9.8 C13.wet 5 13.198 3.554 0.7010 5 9.845
2798 5@9.8 C13.wet 5 14.034 4.378 0.8635 5 9.845
2799 5@9.8 C13.wet 5 14.871 4.549 0.8972 5 9.845
2800 5@9.8 C13.wet 5 15.250 3.913 0.7719 5 9.845
2801 5@9.8 C13.wet 5 15.267 3.890 0.7674 5 9.845
2802 5@9.8 C13.wet 5 15.741 3.706 0.7311 5 9.845
2803 5@9.8 C13.wet 5 15.913 3.680 0.7260 5 9.845
2804 5@9.8 C13.wet 5 16.737 3.756 0.7409 5 9.845
2805 5@9.8 C13.wet 5 17.004 3.439 0.6784 5 9.845
2806 5@9.8 C13.wet 5 17.017 3.440 0.6785 5 9.845
2807 5@9.8 C13.wet 5 17.564 3.359 0.6626 5 9.845
2808 5@9.8 C13.wet 5 18.112 3.167 0.6247 5 9.845
2809 5@9.8 C13.wet 5 18.374 3.084 0.6083 5 9.845
2810 5@9.8 C13.wet 5 18.948 3.137 0.6188 5 9.845
2811 5@9.8 C13.wet 5 19.442 3.182 0.6276 5 9.845
2812 5@9.8 C13.wet 5 19.540 3.187 0.6286 5 9.845
2813 5@9.8 C13.wet 5 20.311 3.399 0.6704 5 9.845
2814 5@9.8 C13.wet 5 21.081 3.354 0.6616 5 9.845
2828 5@9.8 C13.wet 5 28.088 3.054 0.6024 5 9.845
2829 5@9.8 C13.wet 5 29.042 3.193 0.6298 5 9.845
2830 5@9.8 C13.wet 5 29.995 3.152 0.6217 5 9.845
2831 5@9.8 C13.wet 5 30.889 3.450 0.6805 5 9.845
2832 5@9.8 C13.wet 5 31.782 3.611 0.7122 5 9.845
2833 5@9.8 C13.wet 5 32.693 3.711 0.7320 5 9.845
2834 5@9.8 C13.wet 5 33.603 3.642 0.7183 5 9.845
2835 5@9.8 C13.wet 5 34.514 3.401 0.6709 5 9.845
2836 5@9.8 C13.wet 5 35.215 3.578 0.7058 5 9.845
2837 5@9.8 C13.wet 5 35.915 3.648 0.7195 5 9.845
2838 5@9.8 C13.wet 5 36.196 3.623 0.7146 5 9.845
2839 5@9.8 C13.wet 5 36.600 3.729 0.7356 5 9.845
2840 5@9.8 C13.wet 5 37.211 3.703 0.7303 5 9.845
2841 5@9.8 C13.wet 5 37.821 3.531 0.6965 5 9.845
2842 5@9.8 C13.wet 5 38.619 3.223 0.6358 5 9.845
2843 5@9.8 C13.wet 5 39.028 3.077 0.6069 5 9.845
lowCImarker hiCImarker lowCIpos hiCIpos unique
2783 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2784 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2786 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2787 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2788 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2789 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2790 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2791 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2792 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2793 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2794 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2795 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2796 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2797 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2798 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2799 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2800 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2801 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2802 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2803 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2804 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2805 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2806 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2807 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2808 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2809 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2810 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2811 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2812 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2813 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2814 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2828 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2829 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2830 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2831 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2832 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2833 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2834 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2835 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2836 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2837 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2838 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2839 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2840 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2841 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2842 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
2843 C5_896146 C5_11137211 7.41 37.82 5@9.8_C13.wet
new.int<-new.int[new.int$pos<22,]
new.ints[["5@9.8_C13.wet"]]<-c("5@9.8_C13.wet",min(new.int$pos),max(new.int$pos))
new.ints
$`3@58.5_C13.dry`
[1] "3@58.5_C13.dry" "49.2826666666667" "67.6726"
$`5@9.8_C13.wet`
[1] "5@9.8_C13.wet" "5.785" "21.081"
#"5@69.8_C13.wet"
for.lps.sub<-for.lps[for.lps$unique=="5@69.8_C13.wet",]
new.int<-for.lps.sub[for.lps.sub$st.lod>=0.6,]
new.int
qtl.id phenotype chr pos lod st.lod chromosome position
2610 5@69.8 C13.wet 5 68.85 3.645 0.8229 5 69.82
2611 5@69.8 C13.wet 5 69.82 4.429 1.0000 5 69.82
2612 5@69.8 C13.wet 5 70.43 4.414 0.9967 5 69.82
2613 5@69.8 C13.wet 5 71.04 4.230 0.9551 5 69.82
2614 5@69.8 C13.wet 5 71.68 4.022 0.9081 5 69.82
2615 5@69.8 C13.wet 5 72.37 4.064 0.9177 5 69.82
2616 5@69.8 C13.wet 5 73.06 3.914 0.8838 5 69.82
2617 5@69.8 C13.wet 5 73.75 3.594 0.8116 5 69.82
2618 5@69.8 C13.wet 5 73.84 3.451 0.7792 5 69.82
2619 5@69.8 C13.wet 5 74.28 3.490 0.7880 5 69.82
2620 5@69.8 C13.wet 5 74.83 3.458 0.7808 5 69.82
2621 5@69.8 C13.wet 5 75.37 3.291 0.7430 5 69.82
2622 5@69.8 C13.wet 5 76.15 3.314 0.7482 5 69.82
2623 5@69.8 C13.wet 5 76.81 3.250 0.7338 5 69.82
2624 5@69.8 C13.wet 5 77.47 3.090 0.6977 5 69.82
2625 5@69.8 C13.wet 5 78.30 2.740 0.6186 5 69.82
2631 5@69.8 C13.wet 5 81.70 3.119 0.7042 5 69.82
2632 5@69.8 C13.wet 5 82.58 3.695 0.8343 5 69.82
2633 5@69.8 C13.wet 5 83.53 3.518 0.7943 5 69.82
2634 5@69.8 C13.wet 5 84.18 3.608 0.8146 5 69.82
2635 5@69.8 C13.wet 5 84.83 3.235 0.7303 5 69.82
2636 5@69.8 C13.wet 5 85.38 3.256 0.7352 5 69.82
2637 5@69.8 C13.wet 5 85.94 3.245 0.7326 5 69.82
2638 5@69.8 C13.wet 5 86.44 3.405 0.7687 5 69.82
2639 5@69.8 C13.wet 5 87.44 3.429 0.7741 5 69.82
2640 5@69.8 C13.wet 5 88.43 3.397 0.7670 5 69.82
2641 5@69.8 C13.wet 5 89.21 3.790 0.8558 5 69.82
2642 5@69.8 C13.wet 5 89.98 3.982 0.8990 5 69.82
2643 5@69.8 C13.wet 5 90.75 3.926 0.8864 5 69.82
lowCImarker hiCImarker lowCIpos hiCIpos unique
2610 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2611 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2612 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2613 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2614 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2615 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2616 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2617 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2618 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2619 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2620 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2621 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2622 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2623 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2624 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2625 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2631 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2632 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2633 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2634 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2635 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2636 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2637 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2638 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2639 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2640 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2641 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2642 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
2643 C5_21106454 5-26120843 67.88 90.75 5@69.8_C13.wet
new.int<-new.int[new.int$pos<80,]
new.ints[["5@69.8_C13.wet"]]<-c("5@69.8_C13.wet",min(new.int$pos),max(new.int$pos))
new.ints
$`3@58.5_C13.dry`
[1] "3@58.5_C13.dry" "49.2826666666667" "67.6726"
$`5@9.8_C13.wet`
[1] "5@9.8_C13.wet" "5.785" "21.081"
$`5@69.8_C13.wet`
[1] "5@69.8_C13.wet" "68.8515" "78.297"
#"4@2.1_Wilted.dry"
for.lps.sub<-for.lps[for.lps$unique=="4@2.1_Wilted.dry",]
new.int<-for.lps.sub[for.lps.sub$st.lod>=0.6,]
new.int
qtl.id phenotype chr pos lod st.lod chromosome position
1741 4@2.1 Wilted.dry 4 0.0000 2.109 0.7387 4 2.132
1742 4@2.1 Wilted.dry 4 0.7365 2.508 0.8784 4 2.132
1743 4@2.1 Wilted.dry 4 1.4730 2.843 0.9960 4 2.132
1744 4@2.1 Wilted.dry 4 2.1320 2.855 1.0000 4 2.132
1745 4@2.1 Wilted.dry 4 2.7910 2.821 0.9881 4 2.132
1746 4@2.1 Wilted.dry 4 3.5995 2.656 0.9305 4 2.132
1747 4@2.1 Wilted.dry 4 4.4080 2.412 0.8450 4 2.132
1748 4@2.1 Wilted.dry 4 5.3595 2.473 0.8662 4 2.132
1749 4@2.1 Wilted.dry 4 6.3110 2.273 0.7963 4 2.132
1750 4@2.1 Wilted.dry 4 7.2260 1.802 0.6313 4 2.132
1751 4@2.1 Wilted.dry 4 7.8380 1.768 0.6194 4 2.132
1752 4@2.1 Wilted.dry 4 8.0050 1.806 0.6325 4 2.132
1815 4@2.1 Wilted.dry 4 54.6063 1.757 0.6155 4 2.132
1816 4@2.1 Wilted.dry 4 55.3457 1.890 0.6619 4 2.132
1817 4@2.1 Wilted.dry 4 56.0850 1.952 0.6836 4 2.132
1818 4@2.1 Wilted.dry 4 56.4430 2.301 0.8060 4 2.132
1819 4@2.1 Wilted.dry 4 57.2140 2.137 0.7485 4 2.132
1820 4@2.1 Wilted.dry 4 57.7960 2.024 0.7090 4 2.132
1821 4@2.1 Wilted.dry 4 58.3780 1.798 0.6300 4 2.132
1822 4@2.1 Wilted.dry 4 59.2640 2.087 0.7312 4 2.132
1823 4@2.1 Wilted.dry 4 59.7980 2.096 0.7343 4 2.132
1824 4@2.1 Wilted.dry 4 60.3320 1.903 0.6668 4 2.132
1847 4@2.1 Wilted.dry 4 76.0305 1.810 0.6341 4 2.132
1848 4@2.1 Wilted.dry 4 76.8150 1.872 0.6557 4 2.132
1849 4@2.1 Wilted.dry 4 77.1260 1.877 0.6575 4 2.132
1850 4@2.1 Wilted.dry 4 77.4030 1.859 0.6514 4 2.132
1851 4@2.1 Wilted.dry 4 77.7550 1.730 0.6059 4 2.132
lowCImarker hiCImarker lowCIpos hiCIpos unique
1741 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1742 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1743 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1744 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1745 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1746 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1747 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1748 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1749 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1750 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1751 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1752 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1815 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1816 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1817 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1818 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1819 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1820 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1821 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1822 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1823 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1824 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1847 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1848 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1849 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1850 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
1851 4.F6N15 C4_17272166 0 77.75 4@2.1_Wilted.dry
new.int<-new.int[new.int$pos<10,]
new.ints[["4@2.1_Wilted.dry"]]<-c("4@2.1_Wilted.dry",min(new.int$pos),max(new.int$pos))
new.ints
$`3@58.5_C13.dry`
[1] "3@58.5_C13.dry" "49.2826666666667" "67.6726"
$`5@9.8_C13.wet`
[1] "5@9.8_C13.wet" "5.785" "21.081"
$`5@69.8_C13.wet`
[1] "5@69.8_C13.wet" "68.8515" "78.297"
$`4@2.1_Wilted.dry`
[1] "4@2.1_Wilted.dry" "0" "8.005"
#g0
for.lps.sub<-for.lps[for.lps$unique=="1@16.1_g0",]
new.int<-for.lps.sub[for.lps.sub$st.lod>=0.6,]
new.int
qtl.id phenotype chr pos lod st.lod chromosome position
8 1@16.1 g0 1 4.465 3.896 0.6206 1 16.08
23 1@16.1 g0 1 13.136 4.539 0.7229 1 16.08
24 1@16.1 g0 1 13.675 5.180 0.8251 1 16.08
25 1@16.1 g0 1 13.960 5.468 0.8710 1 16.08
26 1@16.1 g0 1 14.832 5.524 0.8798 1 16.08
27 1@16.1 g0 1 15.548 6.049 0.9634 1 16.08
28 1@16.1 g0 1 15.917 6.256 0.9965 1 16.08
29 1@16.1 g0 1 16.076 6.278 1.0000 1 16.08
30 1@16.1 g0 1 16.199 6.251 0.9957 1 16.08
35 1@16.1 g0 1 19.048 3.964 0.6314 1 16.08
36 1@16.1 g0 1 19.710 4.132 0.6582 1 16.08
37 1@16.1 g0 1 19.967 4.444 0.7079 1 16.08
38 1@16.1 g0 1 20.359 4.719 0.7517 1 16.08
39 1@16.1 g0 1 20.927 5.132 0.8174 1 16.08
40 1@16.1 g0 1 21.495 5.321 0.8475 1 16.08
41 1@16.1 g0 1 22.280 5.576 0.8882 1 16.08
42 1@16.1 g0 1 23.064 5.624 0.8959 1 16.08
43 1@16.1 g0 1 23.849 5.459 0.8695 1 16.08
44 1@16.1 g0 1 24.826 4.839 0.7707 1 16.08
45 1@16.1 g0 1 25.380 4.495 0.7160 1 16.08
46 1@16.1 g0 1 25.934 4.102 0.6534 1 16.08
lowCImarker hiCImarker lowCIpos hiCIpos unique
8 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
23 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
24 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
25 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
26 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
27 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
28 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
29 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
30 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
35 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
36 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
37 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
38 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
39 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
40 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
41 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
42 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
43 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
44 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
45 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
46 C1_4617576 C1_7667369 13.96 23.85 1@16.1_g0
new.int<-new.int[new.int$pos>12 & new.int$pos<19.5,]
new.ints[["1@16.1_g0"]]<-c("1@16.1_g0",min(new.int$pos),max(new.int$pos))
new.ints
$`3@58.5_C13.dry`
[1] "3@58.5_C13.dry" "49.2826666666667" "67.6726"
$`5@9.8_C13.wet`
[1] "5@9.8_C13.wet" "5.785" "21.081"
$`5@69.8_C13.wet`
[1] "5@69.8_C13.wet" "68.8515" "78.297"
$`4@2.1_Wilted.dry`
[1] "4@2.1_Wilted.dry" "0" "8.005"
$`1@16.1_g0`
[1] "1@16.1_g0" "13.136" "19.048"
for (i in multi.peaks.unique){
for.rects.sub<-for.rects[for.rects$unique==i,]
for.lps.sub<-for.lps[for.lps$unique==i,]
new.rect<-as.numeric(new.ints[[i]][2:3])
print(ggplot(for.rects.sub, aes(x=pos, y=st.lod, fill=qtl.id, color=qtl.id))+
geom_line(aes(x=pos, y=st.lod),color="black",data=for.lps.sub)+theme_bw()+
geom_rect(aes(xmax=hiCIpos, xmin=lowCIpos, ymin = 0, ymax = 1.02),
color="green",fill=NA)+
geom_rect(xmax=new.rect[2],xmin=new.rect[1],ymin=0,ymax=1.02,
fill="red",color=NA, alpha=0.2)+
ggtitle(paste("lod profile for",i, "\n green box is the original CI, red is the manually set one")))
}
manual.cis <- data.frame(matrix(unlist(new.ints), nrow=6, byrow=T))
manual.cis<-cbind(data.frame(matrix(unlist(strsplit(names(new.ints),"_")), nrow=6, byrow=T)),manual.cis)
manual.cis<-manual.cis[,-3]
manual.cis$chr<-c(3,5,5,4,1,3)
colnames(manual.cis)[1:4]<-c("qtl.id","phenotype","lowCIpos","hiCIpos")
manual.cis<-manual.cis[,c(1,2,5,4,3)]
write.csv(manual.cis, file="dd2014_manualcis.csv", row.names=F)
manual.cis<-read.csv("dd2014_manualcis.csv", header=T)
auto.cis<-newcis.out[,-6]
all.newcis<-rbind(manual.cis,auto.cis)
#combine with untouched cis:
orig.cis<-beststats.ci[-c(1,4,5,6,9,10,12,14,15,18,21),-c(3:5)]
colnames(orig.cis)[2]<-"chr"
all.newcis<-all.newcis[,-1]
all.cis<-rbind(orig.cis,all.newcis)
all.cis<-all.cis[order(all.cis$phenotype,all.cis$chr,all.cis$lowCIpos),]
all.cis<-all.cis[-c(2,7),]
length(all.cis[,1])
[1] 21
write.csv(all.cis, file="dd2014_allcisforcandidates.csv", row.names=F)
toadd<-for.rects[c(4,6),c("phenotype","chr","lowCIpos", "hiCIpos")]
all.cis<-rbind(all.cis,toadd)
all.cis<-all.cis[order(all.cis$phenotype,all.cis$chr,all.cis$lowCIpos),]
colnames(all.cis)<-paste("new",colnames(all.cis), sep="_")
for.rects<-for.rects[order(for.rects$phenotype,for.rects$chr,for.rects$lowCIpos),]
allcistoplot<-cbind(for.rects,all.cis)
ggplot(allcistoplot, aes(x=pos, y=st.lod, fill=qtl.id, color=qtl.id))+
geom_rect(aes(xmax=hiCIpos, xmin=lowCIpos, ymin = 0, ymax = 1.02), color=NA,fill="red", alpha=.2)+
geom_rect(aes(xmax=new_hiCIpos, xmin=new_lowCIpos, ymin = 0, ymax = 1.02), color=NA,fill="blue", alpha=.2)+
geom_line(aes(x=pos, y=st.lod),color="black",data=for.lps)+
facet_grid(phenotype~chromosome, space = "free", scales= "free")+
theme_bw()+
theme(
plot.background = element_blank()
,panel.grid.major = element_blank()
,panel.grid.minor = element_blank(),
strip.background = element_rect(colour="white", fill="white"),
panel.border = element_rect(colour = "black"))+
scale_y_continuous(breaks=c(0,.5,1))+
ggtitle("distribution of standardized LOD profiles overlayed by 90% CIs")
