# To plot a cnv Partition for the PennCNV result, # Bowang Chen # DKFZ # Dec 17, 2010 # if outroot is given, the plot will be saved to png files for each chromosome. plotPennCNV_data<-function(z,outroot=NULL,Chrs=NULL,Colors = NULL, ...){ if (is.null(Colors)) Colors=c("orange3","red","darkblue","green") # colors for opy number =0,1,3,4,2=normal (usually not appear) Samples=as.data.frame(unique(z$SampleID)) colnames(Samples)[1]="SampleID" Samples$index=seq(1:nrow(Samples)) z<-merge(z,Samples,by.x="SampleID",by.y="SampleID") W=8 # inch width H=max(min(nrow(Samples)/50*5,100),5) # height: 30 samples per 5 inches, min 5 inches, max 100 inches. cat("W=",W,", H=",H,"\n") for (i in 1:nrow(z)) { if (z$CN[i]==0) z$Col[i]=Colors[1] else if (z$CN[i]==1) z$Col[i]=Colors[2] else if (z$CN[i]==3) z$Col[i]=Colors[3] else if (z$CN[i]==4) z$Col[i]=Colors[4] else z$Col[i]="black" } if (is.null(Chrs)) { cat(paste("Wrong chromosome set: ",Chrs, "\nSet Chrs=c(1)\n",sep="")) Chrs=c(1) } for (i in Chrs){ cat("plotting chr:",i,", ") if (i>22||i<1) { cat("this is a illegal chr\n") } else { if (is.null(outroot)) { outname=NULL cat("output to screen","\n") } else { if (i<10) outname=paste(outroot,"_chr0",i,".png",sep="") else outname=paste(outroot,"_chr",i,".png",sep="") cat("output to ",outname,"\n") } ch=z[z$Chr==i,] if (nrow(ch)>0){ if (!is.null(outname)) png(outname,units="in",width=W,height=H,res=300) par(xaxt = "s", yaxt = "n",las=0, cex.axis=0.8) plot(ch$StartPosition,ch$index, type = "n", xlab = "",ylab="") rect(ch$StartPosition,ch$index-0.3,ch$EndPosition,ch$index+0.3,border=ch$Col) par(las=2,yaxt="s",cex.axis=0.5) Samples=unique(data.frame(ch$SampleID,ch$index)) axis(2, at = Samples$ch.index, labels = Samples$ch.SampleID) #Y label for sample ID minX=min(as.numeric(ch$StartPosition)) maxX=max(as.numeric(ch$StartPosition)) dx=(maxX-minX)/7 mtext(c(0,1,3,4), side=3, line=1,col=Colors,at=seq(minX+dx*2,minX+dx*5,dx),las=0) # legend of CNV with colors mtext("Sample", 2, line =3, las = 0) mtext(paste("Chromosome ", i,sep=""),1, line = 2.5, las = 0) box() if (!is.null(outname)) dev.off() } } } return(z) } # A raw CNV from PennCNV looks like: # chr1:150822151-150852862 numsnp=32 length=30,712 state2,cn=1 /net/dkfzfsg/gpfs/m/daten/C080-big/Stroke_pennCNV/LRRBAF_new/101_neu startsnp=CN_452248 endsnp=CN_453512 # chr2:89865175-89885829 numsnp=28 length=20,655 state2,cn=1 /net/dkfzfsg/gpfs/m/daten/C080-big/Stroke_pennCNV/LRRBAF_new/101_neu startsnp=CN_867332 endsnp=CN_867355 # chr3:46777921-46788076 numsnp=15 length=10,156 state2,cn=1 /net/dkfzfsg/gpfs/m/daten/C080-big/Stroke_pennCNV/LRRBAF_new/101_neu startsnp=CN_1030794 endsnp=CN_1030808 # chr4:9832514-9843364 numsnp=30 length=10,851 state2,cn=1 /net/dkfzfsg/gpfs/m/daten/C080-big/Stroke_pennCNV/LRRBAF_new/101_neu startsnp=CN_1088266 endsnp=CN_1090406 # chr4:69378740-69408681 numsnp=3 length=29,942 state2,cn=1 /net/dkfzfsg/gpfs/m/daten/C080-big/Stroke_pennCNV/LRRBAF_new/101_neu startsnp=CN_1111686 endsnp=CN_1111690 splitChrIndex<-function(Col=NULL,Sep="=",Index=1){ # to a split raw PennCNV variable into volumns for plot a=strsplit(as.character(Col),Sep) b<-do.call(rbind,a) return=b[,Index] } plotRawPennCNV<-function(filename,outroot=NULL,Chrs=NULL,Colors = NULL, ...){ z<-read.table(filename,header=F) z$Chr=splitChrIndex(z$V1,":",1) z$Chr=splitChrIndex(z$Chr,"r",2) z$Position=splitChrIndex(z$V1,":",2) z$StartPosition=splitChrIndex(z$Position,"-",1) z$EndPosition=splitChrIndex(z$Position,"-",2) z$NumSNP=splitChrIndex(z$V2,"=",2) z$Length=splitChrIndex(z$V3,"=",2) z$State=splitChrIndex(z$V4,",",1) z$State=splitChrIndex(z$State,"e",2) z$CN=splitChrIndex(z$V4,",",2) z$CN=splitChrIndex(z$CN,"=",2) z$SampleID=basename(as.character(z$V5)) z$StartSNP=splitChrIndex(z$V6,"=",2) z$EndSNP=splitChrIndex(z$V7,"=",2) plotPennCNV_data(z,outroot=outroot,Chrs=Chrs,Colors=Colors,...) return(z) } #Windows path format rawname<-"Z:/Stroke_pennCNV/pennCNV_affy_10_neu.adjusted.rawcnv" outroot<-"Z:/Stroke_pennCNV/test" #to plot some chrs, specify them in the Chrs Z<-plotRawPennCNV(rawname,outroot=outroot,Chrs=c(2,3,5)) Z<-plotRawPennCNV(rawname,Chrs=c(2,3)) #linux path format rawname="/net/dkfzfsg/gpfs/m/daten/C080-big/Stroke_pennCNV/pennCNV_affy_10_neu.adjusted.rawcnv" rawname="/net/dkfzfsg/gpfs/m/daten/C080-big/Stroke_pennCNV/LRRBAF_list_76neu_403controls.adjusted.rawcnv" # outroot is the stem of output png files, the plot will be in the form of outroot_chr??.png. outroot<-"/net/dkfzfsg/gpfs/m/daten/C080-big/Stroke_pennCNV/test" # Z<-plotRawPennCNV(rawname,Chrs=c(2,3)) #to plot some chrs, specify them in the Chrs Z<-plotRawPennCNV(rawname,Chrs=c(22,18),outroot=outroot) #to plot all autosomes: Z<-plotRawPennCNV(rawname,Chrs=seq(1,22),outroot=outroot) #to plot some a chr on screen (without png file): Z<-plotRawPennCNV(rawname,Chrs=c(1))