MONOD:LDR2 Calculation: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Shicheng
(Created page with " <nowiki> PairwiseR2plot <- function(R2file=args[1],samplingsize=5000){ file=R2file library(RColorBrewer) require(KernSmooth) data<-read.table(file,sep="\t",skip=1)...")
 
>Shicheng
(Blanked the page)
 
Line 1: Line 1:
<nowiki>
 
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))
</nowiki>

Latest revision as of 19:59, 4 December 2015