R/plotBaseCorrelation.R

In sterlo/kineticWavelets: Small Library for Analysis of PacBio Kinetic Data using wavelets

#'@exportMethod plotBaseCorrelation
setMethod(f='plotBaseCorrelation',"BaseCorrelation",
    function(BaseCorrelation){
    baseCorrelation=BaseCorrelation@baseCorr
    levels=11
    names=colnames(baseCorrelation)
    par(mfcol=c(length(names),length(names)),mar=c(0,0,0,0.5),oma=c(4,3,3,2),xaxt="s",cex=0.5,las=1)
    vals=2^(1:levels)
    for(i in 1:length(names)){
		for(j in 1:length(names)){
		    temp <- vector(length=length(levels));
            if(i==j) {
	            wave1=wd(baseCorrelation[,i],filter.number=1,family="DaubExPhase")
	                for ( level in 1:levels){
	                    temp[level]=sum((accessD(wave1,level=(wave1$nlevels - level)))^2);
	                    }
        		temp <- temp/sum(temp);
     			plot(temp,type="o",col="black",ylim=c(0,0.5),axes=F);
         		box(col="black",lwd=1.5);
	        }
	        if( i < j){
                wave1=wd(baseCorrelation[,i],filter.number=1,family="DaubExPhase")
                wave2=wd(baseCorrelation[,j],filter.number=1,family="DaubExPhase")
	        for ( level in 1:levels){
                test <- cor.test((accessD(wave1,level=(wave1$nlevels - level))),                          		(accessD(wave2,level=(wave2$nlevels - level))));

            temp[level]=test$estimate
            }

            plot(temp,type="l",ylim=c(-1,1),axes=F,xlab="",ylab="");
                   abline(0,0,col="grey",lty=2);
                   box(col="grey"); grid(ny=0,nx=NULL);

            for(level in 1:levels){
                if(temp[level]>0) points(level,temp[level],ty='o',col="red")
                if(temp[level]<0) points(level,temp[level],ty='o',col="blue")
                }

            }

	#smooth top right
            if( i > j){
                wave1=wd(baseCorrelation[,i],filter.number=1,family="DaubExPhase")
                wave2=wd(baseCorrelation[,j],filter.number=1,family="DaubExPhase")
                for ( level in 1:levels){
                    test <- cor.test((accessC(wave1,level=(wave1$nlevels - level))),                          		(accessC(wave2,level=(wave2$nlevels - level))));

                    temp[level]=test$estimate
                }

        plot(temp,type="l",ylim=c(-1,1),axes=F,xlab="",ylab="");
               abline(0,0,col="grey",lty=2);
               box(col="grey"); grid(ny=0,nx=NULL);

                for(level in 1:levels){
                    if(temp[level]>0) points(level,temp[level],ty='o',col="red")
                    if(temp[level]<0) points(level,temp[level],ty='o',col="blue")
                   }
            }
            labels = as.character(2^(1:levels));
            if (i==1) mtext(names[j],side=2,line=0.5,las=3,cex=1.3);
            if (j==1) mtext(names[i],side=3,line=0.5,cex=1.3);
            if (j==length(names)) axis(1,at=1:levels,labels=labels,las=3);
            if (i==length(names) & j==1) axis(4,tck=0.05,cex=0.9) 
            if( i != j ) axis(4,tck=0.05,cex=0.5,labels=FALSE)
            }
    }
});