plugins/base1/se.lu.onk.cgh_dataDumper/trunk/stac_format.R

Code
Comments
Other
Rev Date Author Line
720 11 Jun 08 jari 1 # $Id$
720 11 Jun 08 jari 2
720 11 Jun 08 jari 3 ######################################################################
720 11 Jun 08 jari 4 ################  MAIN PROGRAM #######################################
720 11 Jun 08 jari 5
720 11 Jun 08 jari 6 # Missing values should have the form of Na or NaN
720 11 Jun 08 jari 7 # [reporterId]  [geneSymbol]  [chromosome]  [cytoBand]  [startPosition]  [endPosition]  [assayData = M]
720 11 Jun 08 jari 8
720 11 Jun 08 jari 9 aa<-scan("basematrix.txt",what="character",nlines=1) #aa is header!
720 11 Jun 08 jari 10 cc<-6  #related to format above!
720 11 Jun 08 jari 11 aa.s<-length(aa)-cc #number of samples
720 11 Jun 08 jari 12 my.colClasses<-c("numeric","character","character","character","numeric","numeric",rep("numeric",aa.s))
720 11 Jun 08 jari 13 data<-read.table("basematrix.txt",colClasses=my.colClasses,header=TRUE,na.strings=c(NA))
720 11 Jun 08 jari 14 data$chromosome<-as.integer(data$chromosome)
720 11 Jun 08 jari 15
720 11 Jun 08 jari 16 nbrAssays<-aa.s
720 11 Jun 08 jari 17
720 11 Jun 08 jari 18 assayNames<-read.table("assays.txt",colClasses="character",header=FALSE)
720 11 Jun 08 jari 19
720 11 Jun 08 jari 20 ## sort data and subdivide it
720 11 Jun 08 jari 21 ii<-order(data$chromosome,data$startPosition)
720 11 Jun 08 jari 22 data<-data[ii,] #reorder all data according to chromosom and then to start kb
720 11 Jun 08 jari 23
720 11 Jun 08 jari 24 cloneInfo<-data[,c(seq(1,cc))]
720 11 Jun 08 jari 25 Log2matrix<-as.matrix(data[,-c(seq(1,cc))])
720 11 Jun 08 jari 26 rm(data)
720 11 Jun 08 jari 27   
720 11 Jun 08 jari 28 for(h in 1:length(chromosomes)){ #foreach chromosome
720 11 Jun 08 jari 29   print(h)
720 11 Jun 08 jari 30   workingChrom<-chromosomes[h]
720 11 Jun 08 jari 31   ## Foreach assay in my.assays
720 11 Jun 08 jari 32   for(u in 1:nbrAssays){ #foreach assay
720 11 Jun 08 jari 33     print(assayNames[u])
720 11 Jun 08 jari 34     my.data<-Log2matrix[which(cloneInfo$chromosome==workingChrom),u] #ternary data sorted.
720 11 Jun 08 jari 35     my.info<-cloneInfo[which(cloneInfo$chromosome==workingChrom),]
720 11 Jun 08 jari 36     
720 11 Jun 08 jari 37     ## separate data into p-arm and q-arm, then run each arm independently for gain and losses
720 11 Jun 08 jari 38     p.indexes<-grep("p",my.info$cytoBand)
720 11 Jun 08 jari 39     q.indexes<-grep("q",my.info$cytoBand)
720 11 Jun 08 jari 40     ##
720 11 Jun 08 jari 41     
720 11 Jun 08 jari 42     ## check length of the p.indexes and q.indexes !!! >0 . Bundling not a problem!
720 11 Jun 08 jari 43     if(length(p.indexes)>0){
720 11 Jun 08 jari 44       x1<-c(-99,my.data[p.indexes])
720 11 Jun 08 jari 45       x2<-c(my.data[p.indexes],-99)
720 11 Jun 08 jari 46       z<- x2-x1 #all non-zero values in z corresponds to the breakpoints in the profile
720 11 Jun 08 jari 47       non.zero.z<-which(z !=0) #all breakpoints as c(start1,start2,start3,...,length(data.subset)+1)
720 11 Jun 08 jari 48       
720 11 Jun 08 jari 49       start.i<-non.zero.z[1:length(non.zero.z)-1]
720 11 Jun 08 jari 50       stop.i<-non.zero.z-1
720 11 Jun 08 jari 51       stop.i<-stop.i[-c(1)]
720 11 Jun 08 jari 52       status<-my.data[start.i] #all changes
720 11 Jun 08 jari 53           
720 11 Jun 08 jari 54       ii.g<-which(status==1)
720 11 Jun 08 jari 55       ii.l<-which(status== -1)
720 11 Jun 08 jari 56       
720 11 Jun 08 jari 57       ########
720 11 Jun 08 jari 58       ## write the start stop for the altered states one at a time in a file for each sample
720 11 Jun 08 jari 59     
720 11 Jun 08 jari 60     }
720 11 Jun 08 jari 61     if(length(q.indexes)>0){
720 11 Jun 08 jari 62       x1<-c(-99,my.data[p.indexes])
720 11 Jun 08 jari 63       x2<-c(my.data[p.indexes],-99)
720 11 Jun 08 jari 64       z<- x2-x1 #all non-zero values in z corresponds to the breakpoints in the profile
720 11 Jun 08 jari 65       non.zero.z<-which(z !=0) #all breakpoints as c(start1,start2,start3,...,length(data.subset)+1)
720 11 Jun 08 jari 66       
720 11 Jun 08 jari 67       start.i<-non.zero.z[1:length(non.zero.z)-1]
720 11 Jun 08 jari 68       stop.i<-non.zero.z-1
720 11 Jun 08 jari 69       stop.i<-stop.i[-c(1)]
720 11 Jun 08 jari 70       status<-my.data[start.i] #all changes
720 11 Jun 08 jari 71           
720 11 Jun 08 jari 72       ii.g<-which(status==1)
720 11 Jun 08 jari 73       ii.l<-which(status== -1)    
720 11 Jun 08 jari 74     }
720 11 Jun 08 jari 75     
720 11 Jun 08 jari 76     
720 11 Jun 08 jari 77     #print end resultat to file
720 11 Jun 08 jari 78     
720 11 Jun 08 jari 79   }
720 11 Jun 08 jari 80 }
720 11 Jun 08 jari 81

Generated on Tue Nov 21 03:42:02 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)