MONOD:LDR2 Calculation
Jump to navigation
Jump to search
PairwiseR2plot <- function(R2file=args[1],samplingsize=5000){ file=R2file library(RColorBrewer) require(KernSmooth) data<-read.table(file,sep="\t",skip=1) id<-unique(data[,1]) gap1<-c() gap2<-c() r2<-c() for(j in 1:length(id)){ tmp<-data[data[,1]==id[j],] tmp<-na.omit(tmp) if(nrow(tmp)>2 && sd(tmp[,3],na.rm=T) != 0 && summary(lm(tmp[,3]~tmp[,2],na.action="na.omit"))$coefficients[2,1]<0){ r2<-c(r2,tmp[,3]) gap1<-c(gap1,tmp[,2]) gap2<-c(gap2,tmp[,2]/max(tmp[,2],na.rm=T)) if (length(gap1)>samplingsize) break print(c(j,length(gap1))) } } rlt<-data.frame(gap1,gap2,r2) save(rlt,file=paste(file,".gapR2.RData",sep="")) pdf(paste(file,".paisewise.R2.smooth.pdf",sep="")) par(mfrow=c(2,2)) my.cols <- rev(brewer.pal(11, "RdYlBu")) smoothScatter(x=gap1,y=r2, nrpoints=.3*length(gap1), colramp=colorRampPalette(my.cols), pch=19, cex=.3, col = "green1") smoothScatter(x=gap2,y=r2, nrpoints=.3*length(gap2), colramp=colorRampPalette(my.cols), pch=19, cex=.3, col = "green1") plot(r2~gap1) plot(r2~gap2) dev.off() } file="STL002LG-01.all_chrs.hapInfo.txt.pairwiseR2.gapR2.RData" file="methylC-seq_h1-msc_r2.all_chrs.hapInfo.txt.pairwiseR2.gapR2.RData" file="Colon_primary_tumor.all_chrs.hapInfo.txt.pairwiseR2.gapR2.RData" load(file) gap1<-rlt$gap1 gap2<-rlt$gap2 r2<-rlt$r2 pdf(paste(file,".paisewise.R2.smooth.pdf",sep="")) par(mfrow=c(2,2)) my.cols <- rev(brewer.pal(11, "RdYlBu")) smoothScatter(x=gap1,y=r2, nrpoints=.3*length(gap1), colramp=colorRampPalette(my.cols), pch=19, cex=.3, col = "green1") smoothScatter(x=gap2,y=r2, nrpoints=.3*length(gap2), colramp=colorRampPalette(my.cols), pch=19, cex=.3, col = "green1") plot(r2~gap1) plot(r2~gap2) dev.off() summary(lm(r2~gap1)) summary(lm(r2~gap2))