# TODO: Add comment
# 
# Author: E.Korsching Feb 13, 2012
###############################################################################


acgh.f.score <- function(x, chr.col, phys.pos, threshold=0.1, complete.results=F){
	# according to publication
	# Novel patterns of genome rearrangement and their association with survival in breast cancer
	# James Hicks et al. , Genome Res. 2006 16: 1465-1479 , doi:10.1101/gr.5460106
	# implementation:
	# only those breakpoints with transitions from <threshold to >threshold  or  > (-threshold) to < (-threshold) will be summed up
	# x: dataframe with cols per aCGH experiment, ordered from chr 1 to X/Y, CBS data or similar data structures expected
	# threshold: according to publication but might be adjusted
	# chr.col: chromosome identifier column to x, phys.pos: physical position column to x
	
	#ini
	nr <- nrow(x)
	nc <- ncol(x)
	aCGH.F.score <- matrix(NA, nc, 2)		#contains the number of breakpoints and F scores
	if(complete.results){ bp.postion.mat <- vector("list", nc) }		#contains the breakpoint positions per analysed genome/experiment
	
	#count all breakpoints and segment length in-between
	for(i in 1:nc){
		kk <- 1		#counter 1 to genome length
		ka <- 1		#counter for segment length left to the breakpoint
		ki <- 1		#counter for result matrix
		
		mat.bp <- matrix(NA, nr, 4)	#physical position break point, segment length left, signal value, chromosome
		chr.set <- split(x=x[,i],f=chr.col)		#split experiment into chromosomes
		len.chr.set <- length(chr.set)
		
		for(j in 1:len.chr.set){
			j.len <- length(chr.set[[j]])
			for(k in 1:j.len){
				if(k!=j.len){
					#breakpoint on the first probe position is not counted - if the second signal is also high/low: left segment=1
					#breakpoint on the last probe position of the chr will not be counted, because of missing right segment
					if( (chr.set[[j]][k]<threshold & chr.set[[j]][k+1]>threshold) |					#gain
							(chr.set[[j]][k]>(-threshold) & chr.set[[j]][k+1]<(-threshold)) |		#loss
							(k==1 & chr.set[[j]][k+1]>threshold) |		#we count beginning with the second position, gain
							(k==1 & chr.set[[j]][k+1]<(-threshold))		#loss
					){
						#view: present position to next position
						mat.bp[ki, 1] <- phys.pos[(kk+1)]
						mat.bp[ki, 2] <- ka
						mat.bp[ki, 3] <- x[(kk+1),i]
						mat.bp[ki, 4] <- chr.col[(kk+1)]
						ki <- ki+1
						ka <- 1
#	cat("\n 1 ki ",ki," ka ",ka," kk ",kk)
					}else{
						ka <- ka+1
					}
					kk <- kk+1
#	cat("\n 2 ki ",ki," ka ",ka," kk ",kk)
				}else{	#last position in chr
					mat.bp[ki, 1] <- -1
					if(ka==j.len){
						mat.bp[ki, 2] <- -1		#no breakpoint in chromosome -> filter value below
					}else{
						mat.bp[ki, 2] <- ka		#right segment to last breakpoint
					}
					mat.bp[ki, 3] <- -1		#id chr end, instead signal value
					mat.bp[ki, 4] <- -1
					ki <- ki+1				#stops +1 by last genome position
					ka <- 1
					kk <- kk+1				#stops +1 by last genome/array position
#	cat("\n 3 ki ",ki," ka ",ka," kk ",kk)
				}
			}
		}
		
		#convert result matrix to better format
		if(ki>=3){
			mat.bp <- mat.bp[1:(ki-1), ]	#remove empty matrix elements
			mat.bp <- mat.bp[mat.bp[,2]!=(-1), ]	#remove chromosome(s) without breakpoints
			nrm <- nrow(mat.bp)
			mat.bp1 <- matrix(NA, (nrm-1), 5)	#signal value, left segment length, right segment length, physical position, chromosome
												# -1 because of one liner
			#transform matrix
			ki <- 1		#counter for result matrix
			for(j in 1:(nrm-1)){
				if(mat.bp[j,4]!=(-1)){		#start work only in lines without -1 in col 4
					mat.bp1[j,1] <- mat.bp[j,3]
					mat.bp1[j,2] <- mat.bp[j,2]
					mat.bp1[j,3] <- mat.bp[j+1,2]
					mat.bp1[j,4] <- mat.bp[j,1]
					mat.bp1[j,5] <- mat.bp[j,4]
					ki <- ki+1
				}
			}
			mat.bp1 <- mat.bp1[!is.na(mat.bp1[,5]), ]	#remove empty matrix elements
			nrm <- nrow(mat.bp1)
			if(complete.results){ bp.postion.mat[[i]] <- mat.bp1[ ,c(4,5)] }
			
			#calculate F score of the cases > threshold
			aCGH.F.score[i,1] <- nrm						#number of dected breakpoints
			aCGH.F.score[i,2] <- sum( 2/(mat.bp1[ ,2]+mat.bp1[ ,3]) )			#aCGH F score
			cat("\n number of detected breakpoints obove/below threshold: ", aCGH.F.score[i,1], " aCGH-F-score: ", aCGH.F.score[i,2])
		}else{
			if(complete.results){ bp.postion.mat[[i]] <- 0 }
			cat("\n none detected breakpoints ")
		}
	}
	
	#
	if(complete.results){
		return(list(aCGH.F.score, bp.postion.mat))
	}else{
		return(aCGH.F.score)
	}
}

	
#
#aa <- acgh.f.score(x=HB2.CBS[,7,drop=F], chr.col=HB2.CBS[,1], phys.pos=HB2.CBS[,3], threshold=0.1)
#aa <- acgh.f.score(x=HB2.CBS[,7:72,drop=F], chr.col=HB2.CBS[,1], phys.pos=HB2.CBS[,3], threshold=0.1)

#sum(sapply(aa, length))	# plot as control the breakpoints: marker1