# 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))