## SNPs ## ==== # Asymptotic regression # --------------------- nls.branch1<-nls(boot~SSasympOrig(branch.length,Asym,lrc),data=bootdata,weights=snps.orig) summary(nls.branch1) # bootstrap residuals Eboot<-resid(nls.branch1) # Analysis of bootstrap residuals # ------------------------------- B1<-lm(Eboot~snps.orig*deepness,data=bootdata) summary(B1) plot(B1) # A figure of the regression surface snpseq<-seq(min(bootdata$snps.orig),max(bootdata$snps.orig),len=100) deepseq<-seq(min(bootdata$deepness),max(bootdata$deepness),len=100) colseq<-1-(Eboot-min(Eboot))/(max(Eboot)-min(Eboot)) rsurf.snp<-data.frame(snps.orig=rep(snpseq,100), deepness=rep(deepseq,each=100)) rsurf.snp$pred<-predict(B1,newdata=rsurf.snp) low.colors <- colorRampPalette( c("blue", "lightskyblue1") ) high.colors <- colorRampPalette( c("lightpink1", "red") ) predmat<-matrix(rsurf.snp$pred,nrow=100,ncol=100,byrow=FALSE) predmat.neg<-predmat predmat.neg[predmat.neg>=0]<-NA predmat.pos<-predmat predmat.pos[predmat.pos<0]<-NA par(mar=c(5,5,2,2)) image(deepseq,snpseq,t(predmat.pos),col=high.colors(50), xlab="Node deepness",ylab="Number of SNPs",cex.lab=1.5,cex.axis=1.3) image(deepseq,snpseq,t(predmat.neg),col=low.colors(50),add=TRUE) box(bty="o") contour(deepseq,snpseq,t(predmat),labcex=1.2,method="edge",add=TRUE) points(bootdata$snps.orig~bootdata$deepness,pch=21,bg=gray(colseq),cex=2) ########################################################################### ## Parsimony informative SNPs (PIS) ## ================================ # Asymptotic regression # --------------------- # NOTE! This is the same as above nls.branch1<-nls(boot~SSasympOrig(branch.length,Asym,lrc),data=bootdata,weights=snps.orig) summary(nls.branch1) # bootstrap residuals Eboot<-resid(nls.branch1) # Analysis of bootstrap residuals # ------------------------------- B1.pis<-lm(Eboot~pis.orig*deepness,data=bootdata) summary(B1.pis) plot(B1.pis) # A figure of the regression surface pisseq<-seq(min(bootdata$pis.orig),max(bootdata$pis.orig),len=100) deepseq<-seq(min(bootdata$deepness),max(bootdata$deepness),len=100) colseq<-1-(Eboot-min(Eboot))/(max(Eboot)-min(Eboot)) rsurf.pis<-data.frame(pis.orig=rep(pisseq,100), deepness=rep(deepseq,each=100)) rsurf.bootdata$pred<-predict(B1b,newdata=rsurf.pis) low.colors <- colorRampPalette( c("blue", "lightskyblue1") ) high.colors <- colorRampPalette( c("lightpink1", "red") ) predmat.pis<-matrix(rsurf.bootdata$pred,nrow=100,ncol=100,byrow=FALSE) predmat.pis.neg<-predmat.pis predmat.pis.neg[predmat.pis.neg>=0]<-NA predmat.pis.pos<-predmat.pis predmat.pis.pos[predmat.pis.pos<0]<-NA par(mar=c(5,5,2,2)) image(deepseq,pisseq,t(predmat.pis.pos),col=high.colors(50), xlab="Node deepness",ylab="Number of parsimony informative SNPs",cex.lab=1.5,cex.axis=1.3) image(deepseq,pisseq,t(predmat.pis.neg),col=low.colors(50),add=TRUE) box(bty="o") contour(deepseq,pisseq,t(predmat.pis),labcex=1.2,method="edge",add=TRUE) points(bootdata$pis.orig~bootdata$deepness,pch=21,bg=gray(colseq),cex=2) ########################################################################### ## Loci ## ==== # Asymptotic regression # --------------------- nls.branch2<-nls(boot~SSasympOrig(branch.length,Asym,lrc),data=bootdata,weights=sloci.orig) summary(nls.branch2) # bootstrap residuals Eboot2<-resid(nls.branch3) # Analysis of bootstrap residuals # ------------------------------- B1.loci<-lm(Eboot2~sloci.orig*deepness,data=bootdata) summary(B1.loci) plot(B1.loci) # A figure of the regression surface lociseq<-seq(min(bootdata$sloci.orig),max(bootdata$sloci.orig),len=100) deepseq<-seq(min(bootdata$deepness),max(bootdata$deepness),len=100) rsurf.loci<-data.frame(sloci.orig=rep(lociseq,100), deepness=rep(deepseq,each=100)) rsurf.loci$pred<-predict(B1.loci,newdata=rsurf.loci) colseq2<-1-(Eboot2-min(Eboot2))/(max(Eboot2)-min(Eboot2)) low.colors <- colorRampPalette( c("blue", "lightskyblue1") ) high.colors <- colorRampPalette( c("lightpink1", "red") ) predmat2L<-matrix(rsurf.loci$pred,nrow=100,ncol=100,byrow=FALSE) predmat2L.neg<-predmat2L predmat2L.neg[predmat2L.neg>=0]<-NA predmat2L.pos<-predmat2L predmat2L.pos[predmat2L.pos<0]<-NA par(mar=c(5,5,2,2)) image(deepseq,lociseq,t(predmat2L.pos),col=high.colors(50), xlab="Node deepness",ylab="Number of loci",cex.lab=1.5,cex.axis=1.2) image(deepseq,lociseq,t(predmat2L.neg),col=low.colors(50),add=TRUE) box(bty="o") contour(deepseq,lociseq,t(predmat2L),labcex=1.2,method="edge",add=TRUE) points(bootdata$sloci.orig~bootdata$deepness,pch=21,bg=gray(colseq2),cex=2)