# $Id$

######################################################################
################  MAIN PROGRAM #######################################

# Missing values should have the form of Na or NaN
# [reporterId]	[geneSymbol]	[chromosome]	[cytoBand]	[startPosition]	[endPosition]	[assayData = M]

aa<-scan("basematrix.txt",what="character",nlines=1) #aa is header!
cc<-6  #related to format above!
aa.s<-length(aa)-cc #number of samples
my.colClasses<-c("numeric","character","character","character","numeric","numeric",rep("numeric",aa.s))
data<-read.table("basematrix.txt",colClasses=my.colClasses,header=TRUE,na.strings=c(NA))
data$chromosome<-as.integer(data$chromosome)

nbrAssays<-aa.s

assayNames<-read.table("assays.txt",colClasses="character",header=FALSE)

## sort data and subdivide it
ii<-order(data$chromosome,data$startPosition)
data<-data[ii,] #reorder all data according to chromosom and then to start kb

cloneInfo<-data[,c(seq(1,cc))]
Log2matrix<-as.matrix(data[,-c(seq(1,cc))])
rm(data)
	
for(h in 1:length(chromosomes)){ #foreach chromosome
	print(h)
	workingChrom<-chromosomes[h]
	## Foreach assay in my.assays
	for(u in 1:nbrAssays){ #foreach assay
		print(assayNames[u])
		my.data<-Log2matrix[which(cloneInfo$chromosome==workingChrom),u] #ternary data sorted.
		my.info<-cloneInfo[which(cloneInfo$chromosome==workingChrom),]
		
		## separate data into p-arm and q-arm, then run each arm independently for gain and losses
		p.indexes<-grep("p",my.info$cytoBand)
		q.indexes<-grep("q",my.info$cytoBand)
		##
		
		## check length of the p.indexes and q.indexes !!! >0 . Bundling not a problem!
		if(length(p.indexes)>0){
			x1<-c(-99,my.data[p.indexes])
			x2<-c(my.data[p.indexes],-99)
			z<- x2-x1 #all non-zero values in z corresponds to the breakpoints in the profile
			non.zero.z<-which(z !=0) #all breakpoints as c(start1,start2,start3,...,length(data.subset)+1)
			
			start.i<-non.zero.z[1:length(non.zero.z)-1]
			stop.i<-non.zero.z-1
			stop.i<-stop.i[-c(1)]
			status<-my.data[start.i] #all changes
					
			ii.g<-which(status==1)
			ii.l<-which(status== -1)
			
			########
			## write the start stop for the altered states one at a time in a file for each sample
		
		}
		if(length(q.indexes)>0){
			x1<-c(-99,my.data[p.indexes])
			x2<-c(my.data[p.indexes],-99)
			z<- x2-x1 #all non-zero values in z corresponds to the breakpoints in the profile
			non.zero.z<-which(z !=0) #all breakpoints as c(start1,start2,start3,...,length(data.subset)+1)
			
			start.i<-non.zero.z[1:length(non.zero.z)-1]
			stop.i<-non.zero.z-1
			stop.i<-stop.i[-c(1)]
			status<-my.data[start.i] #all changes
					
			ii.g<-which(status==1)
			ii.l<-which(status== -1)		
		}
		
		
		#print end resultat to file
		
	}
}