| Previous changeset 2:674d50b9ff42 (2016-01-04) Next changeset 4:26bbb1d3ee8a (2016-01-04) |
|
Commit message:
Uploaded |
|
added:
geneDist.R |
| b |
| diff -r 674d50b9ff42 -r 0fc8a9cce0a1 geneDist.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/geneDist.R Mon Jan 04 11:38:53 2016 -0500 |
| [ |
| b'@@ -0,0 +1,394 @@\n+args <- commandArgs()\n+input <- args[4]\n+pngFile <- args[5]\n+dataTable <-read.table(file=input, header=TRUE);\n+chip.data<-data.frame(dataTable)\n+ifReg <- 0\n+if (length(unique(chip.data$Reg))>1) {\n+ ifReg <- 1\n+}\n+\n+ifPDF <- 0\n+bootstrap <- 1\n+\n+if (length(args)>=8) {\n+\tifPDF=args[8]\n+\tbootstrap=args[9]\n+}\n+if (length(args)==7 & args[7]==1) {\n+\tifPDF=1\n+}\n+\n+ifControl <- 0\n+if (length(args)>=7 & args[7]!=1 & args[7]!=0) {\n+ dataTable <-read.table(file=args[7], header=TRUE);\n+ control.data<-data.frame(dataTable)\n+ ifControl <- 1\n+\n+\n+}\n+\n+error.bar <- function(x, y, upper, lower=upper, length=0.1,...){\n+ if(length(x) != length(y) | length(y) !=length(lower) | length(lower) != length(upper))\n+ stop("vectors must be same length")\n+ arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)\n+}\n+\n+\n+\n+logFile <- args[6]\n+sink(logFile, append=FALSE, split=FALSE)\n+\n+if (ifReg & ifControl) {\n+\n+ if (ifPDF==1) {\n+ pdf(file = pngFile, width = 14, height = 13, pointsize = 20, bg = "white")\n+ } else {\n+ png(filename = pngFile, width = 1140, height = 840, units = "px", pointsize = 20, bg = "white", res = NA)\n+ plot(1:10)\n+ }\n+ op <- par(mfrow = c(3,1))\n+} else {\n+\tif (ifReg | ifControl) {\n+\n+\t if (ifPDF==1) {\n+\t\tpdf(file = pngFile, width = 20, height = 17, pointsize = 20, bg = "white")\n+\t } else {\n+\t\tpng(filename = pngFile, width = 1580, height = 1230, units = "px", pointsize = 20, bg = "white", res = NA)\n+\t\tplot(1:10)\n+\t }\n+\t op <- par(mfrow = c(2,1))\n+\t} else {\n+\t if (ifPDF==1) {\n+\t\tpdf(file = pngFile, width = 22, height = 8, pointsize = 20, bg = "white")\n+\t } else {\n+\t\tpng(filename = pngFile, width = 1580, height = 530, units = "px", pointsize = 20, bg = "white", res = NA)\n+\t\tplot(1:10)\n+\t }\n+\t op <- par(mfrow = c(1,1))\n+\t}\n+}\n+myColor <- 1\n+myColor[1] <- colors()[131]\n+myColor[2] <- colors()[59]\n+myColor[3] <- colors()[76]\n+myColor[4] <- colors()[88]\n+myColor[5] <- colors()[17]\n+myColor[6] <- colors()[565]\n+myColor[7] <- colors()[454]\n+myColor[8] <- colors()[401]\n+myColor[9] <- colors()[99]\n+myColorControl <- 1\n+myColorControl[1] <- colors()[24]\n+myColorControl[2] <- colors()[278]\n+myColorControl[3] <- colors()[305]\n+myColorControl[4] <- colors()[219]\n+myColorControl[5] <- colors()[343]\n+myColorControl[6] <- colors()[245]\n+myLevels <- 0\n+nn <- colnames(chip.data)\n+ifk27 <- 0\n+if (length(colnames(chip.data))==41) {\n+\tcc <- c(1:41)\n+\tcolnames(chip.data) <- cc\n+}\n+\n+if (length(colnames(chip.data))==45) {\n+\tcc <- c(1:45)\n+\tcolnames(chip.data) <- cc\n+\tifk27 <- 1\n+}\n+\n+if (!ifk27) {\n+\n+\tcountTotal <- length(unique(chip.data$"1"))\n+\ttt <- which(chip.data$"16">0)\n+\tcountTotalEnh <- length(unique(chip.data$"1"[tt]))\n+\ttt <- which(chip.data$"10">0)\n+\tcountTotalProm <- length(unique(chip.data$"1"[tt]))\n+\ttt <- which(chip.data$"12">0)\n+\tcountTotalImDown <- length(unique(chip.data$"1"[tt]))\n+\ttt <- which(chip.data$"18">0)\n+\tcountTotalIntra <- length(unique(chip.data$"1"[tt]))\n+\ttt <- which(chip.data$"20">0)\n+\tcountTotal5kbD <- length(unique(chip.data$"1"[tt]))\n+\ttt <- which(chip.data$"28">0)\n+\tcountTotal1IntronI <- length(unique(chip.data$"1"[tt]))\n+\ttt <- which(chip.data$"32">0)\n+\tcountTotalExonsI <- length(unique(chip.data$"1"[tt]))\n+\ttt <- which(chip.data$"36">0)\n+\tcountTotalIntronsI <- length(unique(chip.data$"1"[tt]))\n+\ttt <- which(chip.data$"40">0)\n+\tcountTotalJonctionsI <- length(unique(chip.data$"1"[tt]))\n+\n+\tnames <- c("Enh.","Prom.","Imm.Down.","Intrag.","GeneDown.","F.Intron","Exons","2,3,etc.Introns","E.I.Junctions")\n+\tyChIPTotal <- c (countTotalEnh/countTotal,countTotalProm/countTotal, countTotalImDown/countTotal,countTotalIntra/countTotal,countTotal5kbD/countTotal,countTotal1IntronI/countTotal,countTotalExonsI/countTotal,countTotalIntronsI/countTotal,countTotalJonctionsI/countTotal)\n+\n+} else {\n+\tcountTotal <- length(unique(chip.data$"1"))\n+\ttt <- which(chip.data$"42">0)\n+\tcountTSS <- length(unique(chip.data$"1"[tt]))\n+\ttt <- which(chip.data$"44">0)\n+\tcountG'..b'ISubset <- length(unique(subset.data$"1"[tt]))\n+\t\t\t ySubsetTotal <- c (countTotalEnhSubset/countTotalSubset,countTotalPromSubset/countTotalSubset, countTotalImDownSubset/countTotalSubset,countTotalIntraSubset/countTotalSubset,countTotal5kbDSubset/countTotalSubset,countTotal1IntronISubset/countTotalSubset,countTotalExonsISubset/countTotalSubset,countTotalIntronsISubset/countTotalSubset,countTotalJonctionsISubset/countTotalSubset)\n+\t\t\t} else {\n+\t\t\t\ttt <- which(subset.data$"42">0)\n+\t\t\t\tcountTotalTSSSubset <- length(unique(subset.data$"1"[tt]))\n+\t\t\t\ttt <- which(subset.data$"44">0)\n+\t\t\t\tcountTotalGeneSubset <- length(unique(subset.data$"1"[tt]))\n+\t\t\t\tySubsetTotal <- c (countTotalTSSSubset/countTotalSubset,countTotalGeneSubset/countTotalSubset)\n+\t\t\t}\t \n+\t\t\tfor (i in c(1:nlev)) {\t\n+\t\t\tcum[r,i] <- ySubsetTotal[i] \n+\t\t\tcumEnrichTotal[r,i] <- ySubsetTotal[i]/yChIPTotal[i] \n+\t\t\tif (ifControl) {\n+\t\t\t\tcumEnrichControl[r,i] <- ySubsetTotal[i]/yControlTotal[i] \n+\t\t\t} \n+\t\t }\n+\t\t colReg[r]<-myColor[r+3]\n+\t}\n+\tif (ifControl) {\n+\t \tfor (i in c(1:nlev)) {\t\n+\t\t\tcum[4,i] <- yControlTotal[i] \n+\t\t}\t\t\n+\t}\n+\t\n+\tcat ("Proportion of genes with a peak from the ChIP dataset in a given genomic region:\\n")\n+\tfor (r in c(1:length(myLevels))) {\n+\t\tcat (paste(myLevels[r],":\\n",sep=""))\n+\t\tcat (paste(c(names),sep=\'\\t\'))\n+\t\tcat("\\n")\n+\t\tcat (paste(c(cum[r,] ),sep=\'\\t\'))\n+\t\tcat("\\n")\n+\t}\n+\tcat("\\n")\n+\tpar(mar=c(5.1, 7.1, 4.1, 2.1)) \t\n+\tif (ifControl) {\n+\t\tbarx <- barplot(cum,xlab="",beside=TRUE, col=c(colReg, colors()[334]), legend = c(myLevels,"Control"), names.arg=c(names),ylab="Proportion of genes with a peak")\n+\t\tcat ("Proportion of genes with a peak from the Control dataset in a given genomic region:\\n")\n+\t\tcat (paste(c(names),sep=\'\\t\'))\t\n+\t\tcat("\\n")\n+\t\tcat (paste(c(yControlTotal),sep=\'\\t\'))\n+\t\tcat("\\n\\n")\t\t\t\t\n+\t\tif (bootstrap>1) {\n+\t\t\twiskers <- cum-cum\n+\t\t\twiskers[nrow(wiskers),] <- sqrt(colVars)*1.96\n+\t\t\terror.bar(barx,cum,wiskers)\t\n+\t\t\tcat ("Standard deviation for the Control dataset in a given genomic region:\\n")\n+\t\t\tcat (paste(c(names),sep=\'\\t\'))\n+\t\t\tcat("\\n")\n+\t\t\tcat (paste(c(sqrt(colVars)),sep=\'\\t\'))\n+\t\t\tcat("\\n\\n")\t\n+\t\t} \n+\t} else {\n+\t\tbarplot(cum,xlab="",beside=TRUE, col=c(colReg), legend = c(myLevels), names.arg=c(names),ylab="Proportion of genes with a peak")\n+\t}\t\t\n+\tbarplot(cumEnrichTotal-1,xlab="",beside=TRUE, col=c(colReg), names.arg=c(names),ylab="Enrichment in comparison\\nwith the total gene set", yaxt="n")\n+\tminX <- min(cumEnrichTotal-1)\n+\tmaxX <- max(cumEnrichTotal-1)\n+\tx = seq(length=11, from=round(minX*10)/10, by=round((maxX-minX)*10)/100)\n+\taxis(2, at=x,labels=x+1, las=2)\n+\n+\tcat ("Enrichment of genomic regions, Transcriptional categories vs All Genes:\\n")\n+\tfor (r in c(1:length(myLevels))) {\n+\t\t\tcat (paste(myLevels[r],":\\n",sep=""))\n+\t\t\tcat (paste(c(names),sep=\'\\t\'))\t\n+\t\t\tcat("\\n")\n+\t\t\tcat (paste(c(cumEnrichTotal[r,]),sep=\'\\t\'))\t\t\t\n+\t\t\tcat("\\n")\n+\t}\t\n+\tcat("\\n")\n+\n+\tif (ifControl) {\n+\t\tbarplot(cumEnrichControl-1,xlab="",beside=TRUE, col=c(colReg), names.arg=c(names),ylab="Enrichment in comparison\\nwith control", yaxt="n")\n+\t\tminX <- min(cumEnrichControl-1)\n+\t\tmaxX <- max(cumEnrichControl-1)\n+\t\tx = seq(length=11, from=round(minX*10)/10, by=round((maxX-minX)*10)/100)\n+\t\taxis(2, at=x,labels=x+1, las=2)\n+\t\tcat ("Enrichment of genomic regions, ChIP peaks vs Control Peaks:\\n")\n+\t\tfor (r in c(1:length(myLevels))) {\n+\t\t\tcat (paste(myLevels[r],":\\n",sep=""))\n+\t\t\tcat (paste(c(names),sep=\'\\t\'))\t\n+\t\t\tcat("\\n")\n+\t\t\tcat (paste(c(cumEnrichControl[r,]),sep=\'\\t\'))\t\t\t\n+\t\t\tcat("\\n")\n+\t\t}\t\n+\t\tcat("\\n")\n+\t\tif (bootstrap>1) {\n+\t\t\tcat ("Two-side P-values for each genomic category:\\n")\n+\t\t\tfor (r in c(1:length(myLevels))) {\n+\t\t\t\tz <- (cum[r,]-yControlTotal)/sqrt(colVars)\n+\t\t\t\tpvalues <- 2*pnorm(-abs(z))\n+\t\t\t\tcat (paste(myLevels[r],":\\n",sep=""))\n+\t\t\t\tcat (paste(c(names),sep=\'\\t\'))\t\n+\t\t\t\tcat("\\n")\n+\t\t\t\tcat (paste(c(pvalues),sep=\'\\t\'))\t\t\t\n+\t\t\t\tcat("\\n")\n+\t\t\t}\t\n+\t\t}\n+\t}\t\n+} \n+sink() #stop sinking :)\n+dev.off()\n+\n' |