[Patchwork-commits] r170 - in pkg/TAPS: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 3 14:42:06 CEST 2013
Author: sebastian_d
Date: 2013-06-03 14:42:06 +0200 (Mon, 03 Jun 2013)
New Revision: 170
Modified:
pkg/TAPS/DESCRIPTION
pkg/TAPS/R/TAPS.r
pkg/TAPS/man/TAPS_call.Rd
pkg/TAPS/man/TAPS_plot.Rd
Log:
latest TAPS version update to hopefully help R-forge make it build correctly
Modified: pkg/TAPS/DESCRIPTION
===================================================================
--- pkg/TAPS/DESCRIPTION 2013-04-23 14:51:43 UTC (rev 169)
+++ pkg/TAPS/DESCRIPTION 2013-06-03 12:42:06 UTC (rev 170)
@@ -2,9 +2,10 @@
Type: Package
Title: Tumor Abberation Prediction Suite
Version: 1.0
-Date: 2013-02-11
-Author: Markus Rasmussen, Hanna Goransson-Kultima
-Maintainer: Markus Rasmussen <Markus.Mayrhofer at medsci.uu.se>
+Date: 2013-05-11
+Author: Markus Mayrhofer, Hanna Goransson-Kultima, Sebastian DiLorenzo
+Maintainer: Sebastian DiLorenzo <Sebastian.DiLorenzo at medsci.uu.se> #Markus Mayrhofer <Markus.Mayrhofer at medsci.uu.se>
Description: Performs a allele-specific copy number analysis of array data.
License: GPL-2
-Depends: R (>= 2.10), DNAcopy, stats, fields, affxparser
+Depends:R (>= 2.10), DNAcopy, fields, affxparser
+#Imports: DNAcopy, fields, affxparser
\ No newline at end of file
Modified: pkg/TAPS/R/TAPS.r
===================================================================
--- pkg/TAPS/R/TAPS.r 2013-04-23 14:51:43 UTC (rev 169)
+++ pkg/TAPS/R/TAPS.r 2013-06-03 12:42:06 UTC (rev 170)
@@ -1,8 +1,11 @@
-
TAPS_plot <- function(directory=NULL,#xlim=c(-1,2),ylim=c(0,1),
bin=400) {
#Automatically check, and if needed install, packages stats and fields
+#Load stats. It should be in all, at least semi-new, R distributions so we dont need to install.package it or
+#pre-install it
+library(stats)
+
#list.of.packages <- c("stats", "fields")
#new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
#if(length(new.packages)) install.packages(new.packages)
@@ -30,6 +33,7 @@
for (i in 1:length(subs)) try( {
setwd(subs[i])
name <- subs[i]
+ #if (length(grep('SampleData.txt',dir()))==0) save.txt(data.frame(Sample=name,cn2='',delta='',loh='',completed.analysis='no'),file='SampleData.txt')
cat(' ..loading', subs[i])
if(length(grep("*.cyhd.cychp",dir()))==1) ##cyhd sample
@@ -103,7 +107,10 @@
alf=alf[!is.na(alf$Value),]
segments <- readSegments() ## segments if available (CBS recommended)
-
+
+ #Remove NA values that some samples give.
+ segments <- segments[!is.nan(segments$Value),]
+
cat(' ..processing')
if (is.null(segments)) { ## segmentation using DNA-copy if needed (must then be installed)
segments <- segment_DNAcopy(Log2)
@@ -145,7 +152,8 @@
###
###
-TAPS_call <- function(directory=NULL,xlim=c(-1,1),ylim=c(0,1),minseg=1,maxCn=12) {
+TAPS_call <- function(directory=NULL,#xlim=c(-1,1),ylim=c(0,1),
+ minseg=1,maxCn=12) {
## TAPS_call outputs the total and minor allele copy numbers of all segments as a text file, and as images for visual confirmation.
## sampleInfo_TAPS.txt must be present in each sample folder. If TAPS_plot could not make a good guess of the Log-R of copy number 2
## and the Log-R difference to a deletion, you must interpret the scatter plots and edit sampleInfo_TAPS.txt.
@@ -162,10 +170,17 @@
setwd(directory)
#subs <- getSubdirs()
- sampleData <- load.txt('SampleData.txt')
+ if (length(grep('SampleData.txt',dir()))==1)
+ {
+ sampleData <- load.txt('SampleData.txt')
+ }
+ else
+ {
+ sampleData <- load.txt('../SampleData.txt')
+ }
subs=as.character(sampleData$Sample)
#save.txt(sampleData,file='sampleData_old.csv')
-
+
if (is.null(subs)) {
subs=thisSubdir()
setwd('..')
@@ -181,6 +196,15 @@
alf <- readAlf(localDir)
segments <- readSegments()
+ #Some samples throw NA values, we simply remove these.
+ Log2=Log2[!is.nan(Log2$Value),]
+ Log2=Log2[!is.na(Log2$Value),]
+
+ alf=alf[!is.nan(alf$Value),]
+ alf=alf[!is.na(alf$Value),]
+
+ segments <- segments[!is.nan(segments$Value),]
+
segments$Value <- segments$Value-median(Log2$Value)
Log2$Value <- Log2$Value-median(Log2$Value)
@@ -197,7 +221,7 @@
regions=allRegions$regions
save(t,regions,file="regions_t.Rdata")
- karyotype_check(regions$Chromosome,regions$Start,regions$End,regions$log2,regions$imba,regions$Cn,regions$mCn,t,ideogram=NULL,name=name,xlim=xlim,ylim=ylim)
+ karyotype_check(regions$Chromosome,regions$Start,regions$End,regions$log2,regions$imba,regions$Cn,regions$mCn,t,ideogram=NULL,name=name)
karyotype_chromsCN(regions$Chromosome,regions$Start,regions$End,regions$log2,
regions$imba,regions$Cn,regions$mCn,ideogram=NULL,
@@ -271,7 +295,7 @@
segment_DNAcopy <- function(Log2) {
## If segmentation is required, DNAcopy is a good choice. Must be installed.
#library(DNAcopy)
- cat('Using DNAcopy to create segments:\n')
+ cat('..Using DNAcopy to create segments:\n')
segs=NULL
chroms=c(as.character(1:22),'X','Y')
chroms=paste('chr',chroms,sep='')
@@ -1082,7 +1106,7 @@
# }
# }
-karyotype_check <- function(chr,start,end,int,ai,Cn,mCn,t,ideogram=NULL,name='',xlim=c(-1.02,1.02),ylim=0:1) {
+karyotype_check <- function(chr,start,end,int,ai,Cn,mCn,t,ideogram=NULL,name='') { #xlim=c(-1.02,1.02),ylim=0:1) {
## TAPS scatter plot of a full sample, used for visual quality control.
png(paste(name,'.karyotype_check.png',sep=''),width=1300,height=1300)
@@ -1127,12 +1151,12 @@
cex=c(size),
cex.lab=2,
mar=c(0.1,0.1,0.1,0.1),
- main = "",
- xlab = name,
- ylab = "",
- col = col,
- xlim = xlim,
- ylim = ylim)
+ main = "",
+ xlab = name,
+ ylab = "",
+ col = col,
+ xlim=c(-1,2),ylim=c(0,1))
+
text(variants_data[,1],variants_data[,2],
labels=variants,
col='black',
@@ -1493,65 +1517,17 @@
# Added TAPS plot functionality
#------------------------------------------------------------
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
-#Some "what is sent" information that I used to determine which objects and columns were needed to call OverviewPlow
-#karyotype_check(regions$Chromosome,regions$Start,regions$End,regions$log2,
-# regions$imba,regions$Cn,regions$mCn,t,ideogram=NULL,name=name,
-# xlim=xlim,ylim=ylim)
-
-
-
-
-#karyotype <- function(chr,start,end,int,ai,ideogram=NULL,name='',xlim=c(-1.02,1.02),ylim=0:1)
-
-#karyotype(regs$chr,regs$start,regs$end,regs$logs,regs$scores,ideogram=NULL,name=name,xlim=xlim,ylim=ylim)
-
-#karyotype What is sent What it is in the function
-# regs$chr, chr
-# regs$start, start
-# regs$end, end
-# regs$logs, int
-# regs$scores, ai
-# ideogram=NULL, NULL
-# name=name, name
-# xlim=c(-1,2), xlim
-# ylim=c(0,1) ylim
-
-#karyotype_chroms <- function(chr,start,end,int,ai,ideogram=NULL,
-# mchr,mpos,mval,schr,spos,sval,
-# name='',xlim=c(-1.02,1.82),ylim=0:1)
-
-
-#karyotype_chroms(regs$chr,regs$start,regs$end,regs$logs,regs$scores,ideogram=NULL,
-# as.character(Log2$Chromosome),Log2$Start,Log2$Value,as.character(alf$Chromosome),alf$Start,alf$Value,
-# name=name,xlim=xlim,ylim=ylim)
-
-
-#karyotype_chroms What is sent What it is in function
-# regs$chr, chr
-# regs$start, start
-# regs$end, end
-# regs$logs, int
-# regs$scores, ai
-# ideogram=NULL, NULL
-# as.character(Log2$Chromosome), mchr
-# Log2$Start, mpos
-# Log2$Value, mval
-# as.character(alf$Chromosome), schr
-# alf$Start, spos
-# alf$Value, sval
-# name=name, name
-# xlim=c(-1,2), xlim
-# ylim=c(0,1) ylim
-
-
-
-#Example call to OverviewPlot:
-#OverviewPlot(regs$chr,regs$start,regs$end,regs$logs,regs$scores,ideogram=NULL,
-# as.character(Log2$Chromosome),Log2$Start,Log2$Value,as.character(alf$Chromosome),alf$Start,alf$Value,
-# name='wip',xlim=c(-1,2),ylim=c(0,1))
-
OverviewPlot <- function(chr,start,end,int,ai,ideogram=NULL,mchr,mpos,mval,schr,spos,sval,name='',xlim=c(-1,2),ylim=c(0,1))
{
ideogram=getIdeogram()
@@ -1583,23 +1559,6 @@
#This details the margins (left,right,bottom,top) of all the screens.
#See split.screen section of http://seananderson.ca/courses/11-multipanel/multipanel.pdf
#for the most complete explanation.
- # split.screen(as.matrix(data.frame(left=c(0.03,(0.97/8)+0.03,2*(0.97/8)+0.03,3*(0.97/8)+0.03,4*(0.97/8)+0.03,5*(0.97/8)+0.03,6*(0.97/8)+0.03,7*(0.97/8)+0.03,
- # 0.03,(0.97/8)+0.03,2*(0.97/8)+0.03,3*(0.97/8)+0.03,4*(0.97/8)+0.03,5*(0.97/8)+0.03,6*(0.97/8)+0.03,7*(0.97/8)+0.03,
- # 0.03,(0.97/8)+0.03,2*(0.97/8)+0.03,3*(0.97/8)+0.03,4*(0.97/8)+0.03,5*(0.97/8)+0.03,6*(0.97/8)+0.03,7*(0.97/8)+0.03),
- # right=c((0.97/8)+0.03,2*(0.97/8)+0.03,3*(0.97/8)+0.03,4*(0.97/8)+0.03,5*(0.97/8)+0.03,6*(0.97/8)+0.03,7*(0.97/8)+0.03,8*(0.97/8)+0.03,
- # (0.97/8)+0.03,2*(0.97/8)+0.03,3*(0.97/8)+0.03,4*(0.97/8)+0.03,5*(0.97/8)+0.03,6*(0.97/8)+0.03,7*(0.97/8)+0.03,8*(0.97/8)+0.03,
- # (0.97/8)+0.03,2*(0.97/8)+0.03,3*(0.97/8)+0.03,4*(0.97/8)+0.03,5*(0.97/8)+0.03,6*(0.97/8)+0.03,7*(0.97/8)+0.03,8*(0.97/8)+0.03),
- # bottom=c((0.8/3)*2+0.1,(0.8/3)*2+0.1,(0.8/3)*2+0.1,(0.8/3)*2+0.1,(0.8/3)*2+0.1,(0.8/3)*2+0.1,(0.8/3)*2+0.1,(0.8/3)*2+0.1,
- # (0.8/3)+0.1,(0.8/3)+0.1,(0.8/3)+0.1,(0.8/3)+0.1,(0.8/3)+0.1,(0.8/3)+0.1,(0.8/3)+0.1,(0.8/3)+0.1,
- # 0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1),
- # top=c(0.9,0.9,0.9,0.9,0.9,0.9,0.9,0.9,
- # (0.8/3)*2+0.1,(0.8/3)*2+0.1,(0.8/3)*2+0.1,(0.8/3)*2+0.1,(0.8/3)*2+0.1,(0.8/3)*2+0.1,(0.8/3)*2+0.1,(0.8/3)*2+0.1,
- # (0.8/3)+0.1,(0.8/3)+0.1,(0.8/3)+0.1,(0.8/3)+0.1,(0.8/3)+0.1,(0.8/3)+0.1,(0.8/3)+0.1,(0.8/3)+0.1))),1) #screen 1 becomes screen 3 - 26
- # split.screen(as.matrix(data.frame(left=c(0.03,0.03,0.03),
- # right=c(1,1,1),
- # bottom=c(0.1,0.45,0.55),
- # top=c(0.45,0.55,1))),2) #screen 2 become screen 27-29
-
split.screen(as.matrix(data.frame(left=c(rep(seq(from=0.05,to=7*(0.9/8)+0.05,by=(0.9/8)),3)),
right=c(rep(seq(from=(0.9/8)+0.05,to=8*(0.9/8)+0.05,by=(0.9/8)),3)),
bottom=c(rep(2*(0.75/3)+0.15,8),rep((0.75/3)+0.15,8),rep(0.15,8)),
@@ -1687,9 +1646,9 @@
}
#If we are one of the chromosomes plotted to the left edge or bottom, add axis
- if(c %in% c(1,9,17)) axis(side=2,cex.axis=0.6,at=seq(from=0,to=0.8,by=0.2))
+ if(c %in% c(1,9,17)) axis(side=2,cex.axis=0.6,tck=-0.05,at=seq(from=0,to=0.8,by=0.2),las=1)
if(c %in% c(seq(from=17,to=24,by=1))) axis(side=1,cex.axis=0.6,at=seq(from=-0.5,to=1,by=0.5))
- if(c %in% c(8,16,24)) axis(side=4,cex.axis=0.6,at=seq(from=0,to=0.8,by=0.2))
+ if(c %in% c(8,16,24)) axis(side=4,cex.axis=0.6,tck=-0.05,at=seq(from=0,to=0.8,by=0.2),las=1)
#Add titles with sample name (name of folder containing sample) and date
if(c==4)
@@ -1729,8 +1688,21 @@
{
this <- ideogram[ideogram$c==c,]
ix <- chr==this$chr
- mix <- mchr==this$chr & mval>=(-2)
+ mix <- mchr==this$chr #& mval>=(-2)
six <- schr==this$chr
+
+ #Predefine ymin ymax and sequence between them
+ ymin=floor(min(int))-0.5
+ if(ymin > -1)
+ {
+ ymin = -1
+ }
+ ymax=ceiling(max(int))+0.5
+ if(ymax < 1)
+ {
+ ymax = 1
+ }
+ seqminmax=seq(ymin,ymax,by=0.5)
if(c>1)
{
@@ -1753,10 +1725,11 @@
}
#Only plot mval higher than -2
- ix=mval>=(-2)
- mpos = mpos[ix]
- mval = mval[ix]
+ #ix=mval>=(-2)
+ #mpos = mpos[ix]
+ #mval = mval[ix]
+ #Remove every third value to make plotting more sparse in this overview. (saves space and time)
ix = seq(3,to=length(mpos),by=3)
mpos = mpos[-ix]
mval = mval [-ix]
@@ -1771,13 +1744,13 @@
col = '#00000003',
xaxt="n",
axes=F,
- ylim = c(-2,1.5),
+ ylim = c(ymin,ymax),
xlim = c(0,max(mpos))
)
#Add axis to the left & right of signal
- axis(side=2,at=seq(from=-2,to=1.5,by=0.5),cex.axis=0.6,pos=0)
- axis(side=4,at=seq(from=-2,to=1.5,by=0.5),cex.axis=0.6,pos=max(mpos))
+ axis(side=2,tck=-0.025,at=seqminmax,cex.axis=0.6,pos=0,las=1)
+ axis(side=4,tck=-0.025,at=seqminmax,cex.axis=0.6,pos=max(mpos),las=1)
mtext("logRatio",side=2,line=-0.8)
#Add colored segment information for each chromosome, red to blue gradient.
@@ -1791,6 +1764,8 @@
col=rep('#000000',sum(ix))
col[pos[ix] < this$mid] <- colors_p(sum(pos[ix] < this$mid))
col[pos[ix] > this$mid] <- colors_q(sum(pos[ix] > this$mid))
+
+ #Add segments
segments(x0=start[ix],x1=end[ix],
y0=int[ix],y1=int[ix],
col=col,
@@ -1865,9 +1840,9 @@
)
#Add axis to the left,right and below of AI. The below axis is the chromosome numbers 1-24.
- axis(side=2,at=seq(from=0,to=1,by=0.2),cex.axis=0.6,pos=0)
+ axis(side=2,tck=-0.04,at=seq(from=0,to=1,by=0.2),cex.axis=0.6,pos=0,las=1)
axis(side=1,at=pre,labels=c(seq(from="1",to="22"),"X","Y"),cex.axis=0.55)
- axis(side=4,at=seq(from=0,to=1,by=0.2),cex.axis=0.6,pos=max(mpos))
+ axis(side=4,tck=-0.04,at=seq(from=0,to=1,by=0.2),cex.axis=0.6,pos=max(mpos),las=1)
mtext("Allele frequency",side=2,line=-0.8)
mtext("Chromosome numbers",side=1,line=1.5,adj=0.4)
@@ -1887,6 +1862,16 @@
}
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+
#Function for generating individual plots for TAPS_call. Gives whole genome, logRatio, cytoband and allele frequency.
karyotype_chroms <- function(chr,start,end,int,ai,ideogram=NULL,mchr,mpos,mval,schr,spos,sval,name='',xlim=c(-1.02,1.82),ylim=0:1)
{
@@ -2013,7 +1998,7 @@
}
#Insert Y and X axis
- axis(side=2,cex.axis=0.6,tck=0.963,col.ticks='#808080',at=seq(from=0,to=1,by=0.2))
+ axis(side=2,cex.axis=0.6,tck=0.963,col.ticks='#808080',at=seq(from=0,to=1,by=0.2),las=1)
axis(side=1,cex.axis=0.6,tck=0.963,col.ticks='#808080',at=seq(from=-1,to=1.5,by=0.1))
#Titles,date/time and axis labels
@@ -2035,8 +2020,21 @@
par(mgp =c(0.5,0.25,0))
#Select the correct chromosome and remove stuff lower than -1
- mix <- mchr==this$chr & mval>(-1)
+ mix <- mchr==this$chr #& mval>(-1)
+ #Predefine ymin ymax and sequence between them
+ ymin=floor(min(int[ix]))-0.5
+ if(ymin > -1)
+ {
+ ymin = -1
+ }
+ ymax=ceiling(max(int[ix]))+0.5
+ if(ymax < 1)
+ {
+ ymax = 1
+ }
+ seqminmax=seq(ymin,ymax,by=0.5)
+
#Create an index of colors relating to the positions and lengths on this chromosome
col=rep('#000000',sum(ix))
col[pos[ix] < this$mid] <- colors_p(sum(pos[ix] < this$mid))
@@ -2052,7 +2050,7 @@
axes=F,
col = '#00000005',
xlim = c(0,this$length),
- ylim = c(-1,1.5))
+ ylim = c(ymin,ymax))
#Add colored segments based on the logRatio data
segments(x0=start[ix],x1=end[ix],
@@ -2061,12 +2059,14 @@
lwd=4)
#Add X and Y axis. The (side=4) axis is just a black line showing where the data ends
- axis(side=2,tck=0.926,col.ticks='#808080',at=seq(from=-1,to=1.5,by=0.5),
- labels=c("-1","-.5","0",".5","1","1.5"),cex.axis=0.6,pos=0)
- axis(side=4,tck=0,pos=this$length,at=seq(from=-1,to=1.5,by=0.5),
- labels=c("-1","-.5","0",".5","1","1.5"),cex.axis=0.6)
+ axis(side=2,tck=0.926,col.ticks='#808080',at=seqminmax, #seq(from=-1,to=1.5,by=0.5),
+ #labels=c("-1","-.5","0",".5","1","1.5"),
+ cex.axis=0.6,pos=0,las=1)
+ axis(side=4,tck=0,pos=this$length,at=seqminmax,
+ #labels=c("-1","-.5","0",".5","1","1.5"),
+ cex.axis=0.6,las=1)
axis(side=1,tck=0.926,col.ticks='#808080',at=seq(5e6,this$length,by=5e6),
- labels=FALSE,cex.axis=0.6,pos=-1)
+ labels=FALSE,cex.axis=0.6,pos=ymin)
#Add Y axis label
mtext("logRatio",side=2,line=-0.8)
@@ -2141,8 +2141,8 @@
ylim = c(0,1))
#Add X and Y axis as well as a line to the rightmost of the data (side=4)
- axis(side=2,tck=0.926,col.ticks='#808080',at=seq(from=0,to=1,by=0.2),cex.axis=0.6,pos=0)
- axis(side=4,tck=0,at=seq(from=0,to=1,by=0.2),cex.axis=0.6,pos=this$length)
+ axis(side=2,tck=0.926,col.ticks='#808080',at=seq(from=0,to=1,by=0.2),cex.axis=0.6,pos=0,las=1)
+ axis(side=4,tck=0,at=seq(from=0,to=1,by=0.2),cex.axis=0.6,pos=this$length,las=1)
axis(side=1,tck=0.925,col.ticks='#808080',at=seq(5e6,this$length,by=5e6),
labels=paste(seq(5,round(this$length/1e6),5),'',sep=''),cex.axis=0.6,pos=0)
@@ -2156,7 +2156,17 @@
}
}
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+
karyotype_chromsCN <- function(chr,start,end,int,ai,Cn,mCn,ideogram=NULL,mchr,mpos,mval,schr,spos,sval,t,name='',xlim=c(-1.02,1.82),ylim=0:1, maxCn=8)
{
@@ -2296,7 +2306,7 @@
}
#Insert Y and X axis
- axis(side=2,cex.axis=0.6,tck=0.963,col.ticks='#808080',at=seq(from=0,to=1,by=0.2))
+ axis(side=2,cex.axis=0.6,tck=0.963,col.ticks='#808080',at=seq(from=0,to=1,by=0.2),las=1)
axis(side=1,cex.axis=0.6,tck=0.963,col.ticks='#808080',at=seq(from=-1,to=1.5,by=0.1))
#Add cn and mcn labels
@@ -2357,7 +2367,7 @@
lwd=5)
#Add Y axis
- axis(side=2,tck=0.926,col.ticks='#808080',at=seq(0,8,by=1),cex.axis=0.6,pos=0)
+ axis(side=2,tck=0.926,col.ticks='#808080',at=seq(0,8,by=1),cex.axis=0.6,pos=0,las=1)
axis(side=4,labels=F,tck=0,pos=this$length)
#Add Y label and date/time title
@@ -2376,8 +2386,21 @@
par(mgp =c(0.5,0.25,0))
#Select the correct chromosome and remove stuff lower than -1
- mix <- mchr==this$chr & mval>(-1)
+ mix <- mchr==this$chr #& mval>(-1)
+ #Predefine ymin ymax and sequence between them
+ ymin=floor(min(int[ix]))-0.5
+ if(ymin > -1)
+ {
+ ymin = -1
+ }
+ ymax=ceiling(max(int[ix]))+0.5
+ if(ymax < 1)
+ {
+ ymax = 1
+ }
+ seqminmax=seq(ymin,ymax,by=0.5)
+
#Create an index of colors relating to the positions and lengths on this chromosome
col=rep('#000000',sum(ix))
col[pos[ix] < this$mid] <- colors_p(sum(pos[ix] < this$mid))
@@ -2393,7 +2416,7 @@
axes=F,
col = '#00000005',
xlim = c(0,this$length),
- ylim = c(-1,1.5))
+ ylim = c(ymin,ymax))
#Add colored segments based on the logRatio data
segments(x0=start[ix],x1=end[ix],
@@ -2402,11 +2425,12 @@
lwd=4)
#Add X and Y axis. The (side=4) axis is just a black line showing where the data ends
- axis(side=2,tck=0.926,col.ticks='#808080',at=seq(from=-1,to=1.5,by=0.5),
- labels=c("-1","-.5","0",".5","1","1.5"),cex.axis=0.6,pos=0)
+ axis(side=2,tck=0.926,col.ticks='#808080',at=seqminmax, #seq(from=-1,to=1.5,by=0.5),
+ #labels=c("-1","-.5","0",".5","1","1.5"),
+ cex.axis=0.6,pos=0,las=1)
axis(side=4,labels=F,tck=0,pos=this$length)
axis(side=1,tck=0.926,col.ticks='#808080',at=seq(5e6,this$length,by=5e6),
- labels=FALSE,cex.axis=0.6,pos=-1)
+ labels=FALSE,cex.axis=0.6,pos=ymin)
#Add Y axis label
mtext("logRatio",side=2,line=0.3)
@@ -2482,7 +2506,7 @@
ylim = c(0,1))
#Add X and Y axis as well as a line to the rightmost of the data (side=4)
- axis(side=2,tck=0.926,col.ticks='#808080',at=seq(from=0,to=1,by=0.2),cex.axis=0.6,pos=0)
+ axis(side=2,tck=0.926,col.ticks='#808080',at=seq(from=0,to=1,by=0.2),cex.axis=0.6,pos=0,las=1)
axis(side=4,labels=F,tck=0,pos=this$length)
axis(side=1,tck=0.925,col.ticks='#808080',at=seq(5e6,this$length,by=5e6),
labels=paste(seq(5,round(this$length/1e6),5),'',sep=''),cex.axis=0.6,pos=0)
@@ -2497,13 +2521,16 @@
}
}
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
+#------------------------------------------------------------------------------------------------------------------------#
-
-
-
-
-
-
#------------------------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------------------------
#TAPS_region() ; a function for zooming in on a specific region and showing the regions genes.
@@ -2511,8 +2538,6 @@
#------------------------------------------------------------------------------------------------------
-
-#Plot function that returns specifics on a region
TAPS_region <- function(directory=NULL,chr,region,hg18=F)
{
if (is.null(directory))
@@ -2526,510 +2551,575 @@
#Set working directory and make that directorys name the samplename. If it is too long, shorten it.
setwd(directory)
- name = thisSubdir()
- if(nchar(name)>12)
- {
- name = substring(name,1,12)
+ subs <- getSubdirs()
+ if (is.null(subs))
+ { ## check samples = subdirectories or a single sample = current directory
+ subs=thisSubdir()
+ setwd('..')
}
- #From TAPS_plot():
- #save(regs,Log2,alf,segments,file="TAPS_plot_output.Rdata")
+ #store chr value in nchr because i made the newbie mistake of having a second variable also
+ #named chr.
+ nchr = chr
- #Check if TAPS_plot_output.Rdata exists. It must for the plot function to run.
- #In this file are the objects needed to plot.
- if(file.exists("TAPS_plot_output.Rdata")==F)
+ #Loop over all subdirectories, should they exist.
+ for (i in 1:length(subs)) try(
{
- cat("TAPS_region could not find a TAPS_plot_output.Rdata file in the sample directory. \n
- You must run TAPS_plot before TAPS_region. \n")
- }
- else
- {
- load("TAPS_plot_output.Rdata")
- }
+ setwd(subs[i])
+ name <- subs[i]
+ cat(' ..plotting', subs[i],"\n")
- #If chr has been given as character
- if(is.character(chr) == T)
- {
- #If chr is longer than 2 characters ("chr1" but not "12")
- #take out the number and convert it to numeric
- if(nchar(chr)>2)
+ if(nchar(name)>12)
{
- c = strsplit(chr,"chr")
- c = as.numeric(c[[1]][2])
+ name = substring(name,1,12)
}
- #chr has been given as character but without "chr" infront. Convert to numeric.
+
+ #From TAPS_plot():
+ #save(regs,Log2,alf,segments,file="TAPS_plot_output.Rdata")
+
+ #Check if TAPS_plot_output.Rdata exists. It must for the plot function to run.
+ #In this file are the objects needed to plot.
+ if(file.exists("TAPS_plot_output.Rdata")==F)
+ {
+ cat("TAPS_region could not find a TAPS_plot_output.Rdata file in the sample directory. \n
+ You must run TAPS_plot before TAPS_region. \n")
+ }
else
{
- c = as.numeric(chr)
+ load("TAPS_plot_output.Rdata")
}
- }
- #else chr has been given as a numeric and can be used directly.
+
+ #If nchr has been given as character
+ if(is.character(nchr) == T)
+ {
+ #If chr is longer than 2 characters ("chr1" but not "12")
+ #take out the number and convert it to numeric
+ if(nchar(nchr)>2)
+ {
+ c = strsplit(nchr,"chr")
+ c = as.numeric(c[[1]][2])
+ }
+ #chr has been given as character but without "chr" infront. Convert to numeric.
+ else
+ {
+ c = as.numeric(nchr)
+ }
+ }
+ #else chr has been given as a numeric and can be used directly.
+ else
+ {
+ c = nchr
+ }
+
+ #Rename parameters of TAPS_plot_output so we can re-use other plots as template.
+ chr =regs$chr
+ start = regs$start
+ end = regs$end
+ int = regs$logs
+ ai = regs$scores
+ ideogram=NULL
+ mchr = as.character(Log2$Chromosome)
+ mpos = Log2$Start
+ mval =Log2$Value
+ schr = as.character(alf$Chromosome)
+ spos = alf$Start
+ sval =alf$Value
+ xlim=c(-1,2)
+ ylim=c(0,1)
+
+ #check if hg18 is true
+ #then load hg18 knownGene list
+ if(hg18==T)
+ {
+ kg = knownGene_hg18
+ }
else
{
- c = chr
+ kg = knownGene
}
- #Rename parameters of TAPS_plot_output so we can re-use other plots as template.
- chr =regs$chr
- start = regs$start
- end = regs$end
- int = regs$logs
- ai = regs$scores
- ideogram=NULL
- mchr = as.character(Log2$Chromosome)
- mpos = Log2$Start
- mval =Log2$Value
- schr = as.character(alf$Chromosome)
- spos = alf$Start
- sval =alf$Value
- xlim=c(-1,2)
- ylim=c(0,1)
- #check if hg18 is true
- #then load hg18 knownGene list
- if(hg18==T)
- {
- kg = knownGene_hg18
- }
- else
- {
- kg = knownGene
- }
+ #
+ #Here is where pretty normal plotting procedure starts
+ #
+ #Get ideogram
+ ideogram=getIdeogram()
- #
- #Here is where pretty normal plotting procedure starts
- #
+ #Set color gradient for p and q arm of chromosome
+ colors_p <- colorRampPalette(c("#6600FF","#9900CC"),space="rgb")
+ colors_q <- colorRampPalette(c("#CC0099","#CC0000"),space="rgb")
- #Get ideogram
- ideogram=getIdeogram()
+ #filter the data (remove NA) and set some standard parameters
+ ai[is.na(ai)]=0
+ aix=ai!=0
+ chr=chr[aix]
+ start=start[aix]
+ end=end[aix]
+ int=int[aix]
+ ai=ai[aix]
+ pos <- (start+end)/2
+ length=end-start
- #Set color gradient for p and q arm of chromosome
- colors_p <- colorRampPalette(c("#6600FF","#9900CC"),space="rgb")
- colors_q <- colorRampPalette(c("#CC0099","#CC0000"),space="rgb")
+ #Set sizes for circles in whole genome plot
+ size=rep(1,length(chr))
+ size[length>2000000]=2
+ size[length>5000000]=3
+ size[length>10000000]=4
- #filter the data (remove NA) and set some standard parameters
- ai[is.na(ai)]=0
- aix=ai!=0
- chr=chr[aix]
- start=start[aix]
- end=end[aix]
- int=int[aix]
- ai=ai[aix]
- pos <- (start+end)/2
- length=end-start
+ #Select the chromosome
+ this <- ideogram[ideogram$c==c,]
- #Set sizes for circles in whole genome plot
- size=rep(1,length(chr))
- size[length>2000000]=2
- size[length>5000000]=3
- size[length>10000000]=4
+ #Extract regions start an stop
+ #region inputted in form start:stop
+ Rstart = min(region)
+ Rend = max(region)
+
+ #Initialize jpeg
+ jpeg(paste(name,'_',this$chr,'_region_',Rstart,"-",Rend,'.jpeg',sep=''),width=11.7,height=8.3,units="in",res=300)
+
+ #split plot into desired formation
+ split.screen(as.matrix(data.frame(left=c(rep(0.05,4),rep(0.47,3)),
+ right=c(0.45,rep(1,6)),
+ bottom=c(0.50,0.05,0.2,0.3,0.50,0.70,0.75),
+ top = c(0.95,0.2,0.3,0.45,0.70,0.75,0.95)))) #screen 1-5
- #Select the chromosome
- this <- ideogram[ideogram$c==c,]
- #Extract regions start an stop
- #region inputted in form start:stop
- Rstart = min(region)
- Rend = max(region)
-
- #Initialize jpeg
- jpeg(paste(name,'_',this$chr,'_region_',Rstart,"-",Rend,'.jpeg',sep=''),width=11.7,height=8.3,units="in",res=300)
-
- #split plot into desired formation
- split.screen(as.matrix(data.frame(left=c(rep(0.05,4),rep(0.47,3)),
- right=c(0.45,rep(1,6)),
- bottom=c(0.50,0.05,0.2,0.3,0.50,0.70,0.75),
- top = c(0.95,0.2,0.3,0.45,0.70,0.75,0.95)))) #screen 1-5
+ #Screen configuration overview
+ #-------------------------------------------------------------
+ #-------------------------------------------------------------
+ # | |
+ # | 7 |
+ # | |
+ # |------------------------------|
+ # 1 | 6 |
+ # |------------------------------|
+ # | |
+ # | 5 |
+ # | |
+ #-----------------------------|------------------------------|
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/patchwork -r 170
More information about the Patchwork-commits
mailing list