Name:
Location: United States

Wednesday, June 21, 2006

DNA pattern: biggest cluster

intervalSizes = seq(from=500, to=12000, by=500)
possibleRepOrigin = data.frame(intervalSize = intervalSizes, begin=rep(NA, length(intervalSizes)), end=rep(NA, length(intervalSizes)), unitClustSize= rep(NA, length(intervalSizes)))

for(i in 1:length(intervalSizes)) {
size = intervalSizes[i]
intervalHead = seq(from=1, to=dnalen, by=size)
if((intervalHead[length(intervalHead)]+size-1)>dnalen) {
numIntervals = length(intervalHead)-1
} else {
numIntervals = length(intervalHead)
}
palCountsEachInterval = rep(NA, numIntervals)
for(intv in 1:numIntervals) {
palCountsEachInterval[intv] = sum(intervalHead[intv]:(intervalHead[intv]+size-1) %in% posvec)
}
possibleRepOrigin$unitClustSize[i] = max(palCountsEachInterval)/size
possibleRepOrigin$begin[i] = intervalHead[which(palCountsEachInterval==max(palCountsEachInterval))[1]] # in case of a tie
if((possibleRepOrigin$begin[i]+size-1)>dnalen) {
possibleRepOrigin$end[i] = dnalen
} else {
possibleRepOrigin$end[i] = possibleRepOrigin$begin[i]+size
}
}
print(possibleRepOrigin)