source: plugins/base1/se.lu.onk.cgh_dataDumper/stac_format.R @ 720

# $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
    
  }
}