# R/find.ladder.R In Fragman: Fragment Analysis in R

```find.ladder <-function(x, ladder, ci.upp=1.96, ci.low=1.96, draw=TRUE, dev=50, warn=TRUE, init.thresh=NULL, sep.index=8, method=NULL, avoid=1000, who="sample", attempt=10, cex.title=0.8){
#roundUP <- function(x) 10^ceiling(log10(x))

hohoho <- big.peaks.col(x[1:length(x)], 100)

if(is.null(method)){
## check if the ladder passes the length of x of is a small ladder
tototo <- big.peaks.col(x[1:length(x)], median(hohoho\$hei)*2)
close <- (length(x) - tototo\$pos[length(tototo\$hei)])
if(close < 100){
method="iter"
}else{
method="iter2"
}
}

if(is.null(init.thresh)){
if(mean(x) > 1000){
#hohoho2 <- big.peaks.col(x[1:length(x)], 2*median(hohoho\$hei))
init.thresh <- quantile(hohoho\$hei,.95)
}else{
init.thresh <- median(hohoho\$hei)/2.2#quantile(hohoho\$hei,.32)
}
}

MSE <- function(x,y){
X <- cbind(1, x)
qr.X <- qr(X)
b <- t(qr.Q(qr.X)) %*% y
R <- qr.R(qr.X)
beta <- as.vector(backsolve(R, b))
fit <- X%*%beta
res <- list(mse=sum((y-fit)^2), beta=beta[2])
return(res)
}
MSE2 <- function(x,y){
X <- cbind(1, x)
qr.X <- qr(X)
b <- t(qr.Q(qr.X)) %*% y
R <- qr.R(qr.X)
beta <- as.vector(backsolve(R, b))
fit <- X%*%beta
mse <- sum((y-fit)^2)# is actually SSE
sst <- sum((y-mean(y))^2)
r2 <- 1 - (mse/sst) #cor(x,y)^2#mse/sum((y-mean(y,na.rm=TRUE))^2)
res <- list(mse=mse, beta=beta, r2=r2)
return(res)
}
MSE3 <- function(mix,miy){
X <- cbind(1, mix)
qr.X <- qr(X)
b <- t(qr.Q(qr.X)) %*% miy
R <- qr.R(qr.X)
beta <- as.vector(backsolve(R, b))
fit <- X%*%beta
mse <- sum((miy-fit)^2)# is actually SSE
sst <- sum((miy-mean(miy))^2)
r2 <- 1 - (mse/sst) #cor(x,y)^2#mse/sum((y-mean(y,na.rm=TRUE))^2)
res <- list(mse=mse, beta=beta, r2=r2)
return(res)
}
##### stablish a minimum threshold
thresh=init.thresh
roxy <- big.peaks.col(x[1:length(x)], thresh) # find all peaks with that threshold

nnn <- length(roxy\$pos)

fdsa=1
#print(paste("Warning: Your ladder for the sample",attr(x,"name") ,"is bad, very low intensity in RFU was found, we are reducing the threshold 2x trying to find it, don't worry too much the analysis will keep going anyways"))

if(fdsa == 1){
cat("\nReducing threshold 2x to find ladder \n")
#cat("There's too many ladder peaks expected and not many peaks found in this sample\nPlease make sure all your ladder peaks specified show up in the sample\nOtherwise provide a smaller ladder so all peaks can be found\n")
}
thresh=thresh/2
roxy <- big.peaks.col(x[1:length(x)], thresh) # find all peaks with that threshold
nnn <- length(roxy\$pos)
fdsa<-fdsa+1
}
####
#use the indicator peak at the beggining which is usually the tallest peak in the ladder
whot <- length(roxy\$pos)*.2
what <- which(roxy\$hei == max(roxy\$hei))
#what <- which(roxy\$hei == sort(roxy\$hei, decreasing = TRUE)[2])
# if you tallest peak is far beyond don't adjust
#if(what > whot){
#  roxy <- roxy
#}else{ max(roxy\$hei)
roxy <- separate(roxy, shift=sep.index, type="pos")

#} plot(x, type="l",xlim=c(2000,7000)); abline(v=roxy\$pos, col="red", lty=3)
#####################
if (method == "iter") {
ii <- which(roxy\$hei == max(roxy\$hei)) + 1
iii <- length(roxy\$hei)
roxy <- list(pos = roxy\$pos[ii:iii], hei = roxy\$hei[ii:iii])

step1 <- combn(roxy\$pos[1:attempt], 3)
step2 <- apply(step1/10, 2, MSE, y = ladder[1:3])
mse <- unlist(lapply(step2, function(x) {
x\$mse
}))
covs <- apply(step1, 2, function(x, y) {
cov(x, y)
step2 <- mse * covs
step3 <- step1[, which(step2 < sort(step2, decreasing = FALSE)[20])]
step4 <- apply(step3, 2, function(x, y) {
which(y %in% x)
}, y = roxy\$pos)
caller <- function(roxy, www, ladder.call, x) {
threshold <- length(x)
posi <- numeric()
expect <- roxy\$pos[www]
modx <- lm(expect ~ poly(xxx, degree = 1))
expecto <- predict(modx, data.frame(xxx = ladder.call))
0.85)]
if (available > 0) {
if ((length(ladder.call) - 1) < 3) {
}
else {
}
}
else {
tope <- ava2 - 2
}
expect <- rep(NA, tope + 1)
for (i in 3:tope) {
if (i == 3 & i != tope) {
expect[1:3] <- roxy\$pos[www]
mod <- MSE2(xxx, expect[1:3])
beta <- (mod)\$beta
expecto <- as.vector(beta[1] + matrix(ladder.call) %*%
beta[-1])
act <- roxy\$pos[-which(roxy\$pos %in% expect)]
yoyo <- abs(expecto[i + 1] - act)
good <- which(yoyo == min(yoyo,na.rm = TRUE))
expect[i + 1] <- act[good]
if (mod\$r2 < 0.9) {
i = tope
}
}
if (i > 3 & i <= 5) {
mod <- MSE2(xx, expect[1:i])
beta <- (mod)\$beta
expecto <- as.vector(beta[1] + matrix(ladder.call) %*%
beta[-1])
act <- roxy\$pos[-which(roxy\$pos %in% expect)]
yoyo <- abs(expecto[i + 1] - act)
#print(min(yoyo))
good <- which(yoyo == min(yoyo, na.rm = TRUE))
expect[i + 1] <- act[good]
if (mod\$r2 < 0.9) {
i = tope
}
}
if (i > 5) {
mod <- MSE2(xx, expect[1:i])
beta <- (mod)\$beta
if (length(which(is.na(beta))) > 0) {
beta[which(is.na(beta))] <- 0
}
expecto <- cbind(rep(1, dim(toto)[1]), toto) %*%
beta
act <- roxy\$pos[-which(roxy\$pos %in% expect)]
yoyo <- abs(expecto[i + 1] - act)
good <- which(yoyo == min(yoyo,na.rm = TRUE))
expect[i + 1] <- act[good]
if (is.na(mod\$r2)) {
mod\$r2 <- 0.1
}
if (mod\$r2 < 0.9) {
i = tope
}
}
if (i == tope & i != 3) {
if (i < 5) {
expect[1:3] <- roxy\$pos[www]
}
else {
}
mod <- MSE2(xx, expect[1:i])
beta <- (mod)\$beta
if (length(which(is.na(beta))) > 0) {
beta[which(is.na(beta))] <- 0
}
if (i < 5) {
}
else {
}
expecto <- cbind(rep(1, dim(toto)[1]), toto) %*%
beta
act <- roxy\$pos[-which(roxy\$pos %in% expect)]
yoyo <- abs(expecto[i + 1] - act)
good <- which(yoyo == min(yoyo, na.rm = TRUE))
expect[i + 1] <- act[good]
}
if (i == tope & i == 3) {
expect[1:3] <- roxy\$pos[www]
}
}
posi <- expect
tutu <- abs(length(x) - posi)
posi <- posi[1:which(tutu == min(tutu, na.rm = TRUE))]
heii <- roxy\$hei[which(roxy\$pos %in% posi)]
fact3 <- length(posi)/fact2
if (length((posi)) < 6) {
poly(posi, degree = length((posi)) - 1)))\$r.squared *
fact3
}
else {
poly(posi, degree = 5)))\$r.squared * fact3
}
roxy2 <- list(pos = posi, hei = heii, wei = ladder.call[1:length(posi)],
error = fact)
return(roxy2)
}
rt <- apply(data.frame(step4), 2, FUN = caller, roxy = roxy,
corrs3 <- unlist(lapply(rt, function(x) {
x\$error
}))
roxy3 <- rt[[which(corrs3 == max(corrs3))]]
if (draw == TRUE) {
limi <- sort(roxy3\$hei, decreasing = TRUE)
plot(x, type = "l", xaxt = "n", ylim = c(0, (limi[3] +
1000)), cex.axis = 0.6, las = 2, xlim = c((min(roxy3\$pos) -
100), (max(roxy3\$pos) + 100)), col = transp("grey35",
0.7), ylab = "RFU", xlab = "", lwd = 2, main = attributes(x)\$mycomm,
cex.main = cex.title)
axis(1, at = roxy3\$pos, labels = roxy3\$wei, cex.axis = 0.6)
points(x = roxy3\$pos, y = roxy3\$hei, cex = 1.1, col = transp("black",
0.85))
points(x = roxy3\$pos, y = roxy3\$hei, pch = 20, col = transp("red",
0.7))
legend("topleft", legend = paste("Correlation:",
round(roxy3\$corr, digits = 4), sep = ""), bty = "n")
legend("topright", legend = c("Peaks selected"),
col = c("red"), bty = "n", pch = c(20), cex = 0.85)
}
roxy <- roxy3
}
if(method == "iter2"){

#caller
last <- length(roxy\$pos)
if((last-attempt) < 0){
roxy\$corr <- 0
roxy\$error <- 0
if(draw == TRUE){
limi <- sort(roxy\$hei, decreasing = TRUE)
plot(x, type="l", xaxt="n", ylim=c(0,(limi[3]+1000)), cex.axis=0.6, las=2,  col=transp("grey35",0.7), ylab="RFU", xlab="", lwd=2, main=attributes(x)\$mycomm, cex.main=cex.title)
axis(1, at=roxy\$pos, labels=roxy\$wei, cex.axis=0.6)
points(x=roxy\$pos, y=roxy\$hei,cex=1.1, col=transp("black",0.85))
points(x=roxy\$pos, y=roxy\$hei, pch=20, col=transp("red",0.7))
legend("topleft", legend=paste("Correlation:",round(roxy\$corr, digits=4), sep=""), bty="n")
legend("topright", legend=c("Peaks selected"), col=c("red"), bty = "n", pch=c(20), cex=0.85)
legend("center", legend=c("Intensity too low to be detected"), col=c("red"), bty = "n", cex=0.85)
}
}else{
# by default attempt=10 making all combinations of 1st 10 peaks
step1 <- combn(roxy\$pos[last:(last-attempt)],3)
#### check MSE
step2 <- apply(step1/10, 2, MSE, y=ladder[lld:(lld-2)]) # function(x,y){sum((summary(lm(I(x)~y))\$residuals)^2)}create a custom function with matrices to extract residuals and calculate MSE
#### mean squared errors
mse <- unlist(lapply(step2, function(x){x\$mse}))
##### covariances
covs <- apply(step1, 2, function(x,y){cov(x,y)}, y=ladder[lld:(lld-2)])
##### scores for selection
step2 <- mse #* covs
##### selected based on scores
step3 <- step1[,which(step2 < sort(step2, decreasing=FALSE)[20])] # 15 models with least MSE
# winner = 1
#for(i in 1:dim(step3)[2]){ # see selected
#  plot(x,type="l", main=i)
#  abline(v=step3[,i], lty=3, col="red")
#}
# position of peaks-combinations selected
step4 <- apply(step3,2,function(x,y){sort(which(y %in% x), decreasing = TRUE)}, y=roxy\$pos) # which peaks are
#caller  6
############
#-----------
############
threshold <- length(x)

posi <- numeric()
############################################################
## short initial model to avoid people using a ladder.call too long, we cut the ladder.call if necessary
expect <- roxy\$pos[www] # initial positions for this combo (last)
modx <- lm(expect~poly(xxx, degree=1))
## assuming model is correct predict where should be the rest of the ladder
expecto <- predict(modx, data.frame(xxx=ladder.call))# + facto

available <- length(roxy\$pos) - length(ladder.call) # there's x extra peaks
ava2 <-  length(ladder.call) - abs(available) # is it more than the user specified??
if(available < 0){ # if there is less peaks available than existing in ladder BAD USER
tope <- 3 # is bad sample

}else{
} # if there is more peaks available than expecting in ladder GOOD USER

#}else{tope <- abs(available)-1}
#####

for(i in 3:(tope-1)){ # for all possible peaks
if(i == 3 & i != tope){# initial step

expect[tope:(tope-2)] <- roxy\$pos[www]
#mod1 <- lm(expect[1:3]~poly(xxx, degree=1))
#expecto <- predict(mod1, data.frame(xxx=ladder.call))# + facto
mod <- MSE3(mix=xxx, miy=expect[tope:(tope-2)])
beta <- (mod)\$beta # extract intercept and slope
## predict the rest of the ladder
expecto <- as.vector(beta[1] + matrix(sort(ladder.call, decreasing = TRUE)) %*% beta[-1])
# remove missing data
condo <- sort(expect[which(!is.na(expect))], decreasing = TRUE)
# keep roxy values that were picked already
act <- sort(roxy\$pos[-which(roxy\$pos %in% condo)], decreasing = TRUE)
# see which is the next most likely value
yoyo <- abs(expecto[i+1] - act)
good <- which(yoyo == min(yoyo, na.rm = TRUE))

qwer <- i
qwer2 <- length(expect)-qwer
expect[qwer2] <- act[good]
# --------------------------------
#if(summary(mod1)\$r.squared < .9){
if(mod\$r2 < .9){
}
# --------------------------------
}else if(i > 3 & i <= 5 ){

#mod1 <- lm(expect[1:i]~poly(xx, degree=1))

mod <- MSE3(mix=xx, miy=expect[tope:qwer2])
beta <- (mod)\$beta # extract intercept and slope
## predict the rest of the ladder
expecto <- as.vector(beta[1] + matrix(sort(ladder.call, decreasing = TRUE)) %*% beta[-1])
# remove missing data
condo <- sort(expect[which(!is.na(expect))], decreasing = TRUE)
# keep roxy values that were picked already
act <- sort(roxy\$pos[-which(roxy\$pos %in% condo)], decreasing = TRUE)
# see which is the next most likely value
yoyo <- abs(expecto[i+1] - act)
good <- which(yoyo == min(yoyo, na.rm = TRUE))
qwer <- i
qwer2 <- length(expect)-qwer
expect[qwer2] <- act[good]
# --------------------------------
#if(summary(mod1)\$r.squared < .9){
if(mod\$r2 < .9){
}
# --------------------------------
}else if(i > 5 ){
#mod1 <- lm(expect[1:i]~poly(xx, degree=4, raw=TRUE))
expect2 <- sort(expect, decreasing = TRUE)
mod <- MSE3(mix=xx, miy=expect2)
beta <- (mod)\$beta
if(length(which(is.na(beta))) > 0){
beta[which(is.na(beta))] <- 0
}
expecto <- cbind(rep(1,dim(toto)[1]),toto) %*% beta

# remove missing data
condo <- sort(expect[which(!is.na(expect))], decreasing = TRUE)
# keep roxy values that were picked already
act <- sort(roxy\$pos[-which(roxy\$pos %in% condo)], decreasing = TRUE)
# see which is the next most likely value
yoyo <- abs(expecto[i+1] - act)
good <- which(yoyo == min(yoyo, na.rm = TRUE))
qwer <- i
qwer2 <- length(expect)-qwer
expect[qwer2] <- act[good]
# --------------------------------
#if(summary(mod1)\$r.squared < .9){
if(is.na(mod\$r2)){ # if model is too bad just assign a zero value
mod\$r2 <- 0.1
}
if(mod\$r2 < .9){
}
# --------------------------------
}else if(i == tope & i == 3){
expect[1:3] <- roxy\$pos[www]
}

##
} # end of for loop
#plot(x, type="l",ylim=c(0,800)); abline(v=expect, col="red", lty=3)
###########################################
posi <- expect
## get rid of selected peaks after reaching the maximum values
#tutu <- abs(length(x) - posi)
#if(length(posi) > 3){
#posi <- posi[1:which(tutu == min(tutu,na.rm = TRUE))]
#}
heii <- roxy\$hei[which(roxy\$pos %in% posi)]

fact3 <- length(posi)/fact2

if(length((posi)) < 6){ # good user
fact <- summary(lm(ladder.call[1:length(posi)]~poly(posi, degree=length((posi))-1)))\$r.squared * fact3
}else{
fact <- summary(lm(ladder.call[1:length(posi)]~poly(posi, degree=5)))\$r.squared * fact3
}

return(roxy2)
}
############
#----------- END OF CALLER
############
## end of caller
# www <- data.frame(step4)[,1]
#for(g in 1:dim(step4)[2]){
#  plot(x,type="l",ylim=c(0,1000))
#}

corrs3 <- unlist(lapply(rt, function(x){x\$error})) #; dis[which(dis == Inf)] <- 1
#corrs <- unlist(lapply(rt, function(x){x[[4]]})) # extract correlations
roxy3 <- rt[[which(corrs3 == max(corrs3))[1]]]
#roxy3 <- rt[[5]]
# plot(x, type="l",xlim=c(2000,7000), ylim=c(0,20000)); abline(v=roxy3\$pos, col="red")
if(draw == TRUE){
limi <- sort(roxy3\$hei, decreasing = TRUE)
plot(x, type="l", xaxt="n", ylim=c(0,(limi[3]+1000)), cex.axis=0.6, las=2,  col=transp("grey35",0.7), ylab="RFU", xlab="", lwd=2, main=attributes(x)\$mycomm, cex.main=cex.title)
axis(1, at=roxy3\$pos, labels=roxy3\$wei, cex.axis=0.6)
points(x=roxy3\$pos, y=roxy3\$hei,cex=1.1, col=transp("black",0.85))
points(x=roxy3\$pos, y=roxy3\$hei, pch=20, col=transp("red",0.7))
legend("topleft", legend=paste("Correlation:",round(roxy3\$corr, digits=4), sep=""), bty="n")
legend("topright", legend=c("Peaks selected"), col=c("red"), bty = "n", pch=c(20), cex=0.85)

}
roxy <- roxy3
}
}
#####################
#####################
#####################
if(method == "cicor"){
roxy <- separate(roxy, shift=40, type="pos")
mm <- median(roxy\$hei)
vvv <- which(roxy\$hei == max(roxy\$hei))
se <- sd(roxy\$hei[-vvv])/sqrt(length(roxy\$hei[-vvv]))
reduced <- roxy\$pos[which(roxy\$hei < mm+(1*se) & roxy\$pos > mm-(1*se) & roxy\$hei > init.thresh)]
reduced2 <- roxy\$pos[which(roxy\$pos >= min(reduced))]

## if there's still just toomany roxy after the reduced search
if(length(reduced2) >= nnn & length(reduced2) < nnn+8){
cors <- apply(all.combs, 2, function(x,y){cor(x,y)}, y=ladder)
# positions of roxy found
found <- all.combs[,which(cors == max(cors))]

nono <- length(x) - max(roxy\$pos)

if(draw == TRUE){
plot(x, type="l", xaxt="n", ylab="Intensity", col=transp("gray29",0.6), main=attributes(x)\$mycomm, cex.main=cex.title)
axis(1, at=roxy\$pos, labels=roxy\$wei, cex.axis=0.6)
abline(v=found, col="red", lty=3)
}
# reduced search to maximum correlations
}else{# there's actually no ladder so no good correlation was found
roxy <- list(pos=seq(1,nnn) + rnorm(nnn,0,1), hei=seq(1,nnn)+ rnorm(nnn,0,1), wei= ladder)
print("Friend I was not able to find a ladder in this sample or you used the wrong ladder, look the plot")
plot(x, type="l", col=transp("black",0.5), main=attributes(x)\$mycomm, cex.main=cex.title)
}
}
#####################
if(method == "red"){

if(nnn > 15){
yoy2 <- rev(seq(1,yoy,by=1))
corres <- numeric()
for(i in 1:yoy){
}
vv <- which(corres == max(corres))[1]
#lapply(roxy,length)
if(draw == TRUE){
plot(x, type="l", xaxt="n", ylim=c(0,(max(roxy\$hei)+1000)), cex.axis=0.6, las=2, xlim=c((min(roxy\$pos)-100),(max(roxy\$pos)+100)), col=transp("grey35",0.7), ylab="RFU", xlab="", lwd=2, main=attributes(x)\$mycomm,cex.main=cex.title)
axis(1, at=roxy\$pos, labels=roxy\$wei, cex.axis=0.6)
points(x=roxy\$pos, y=roxy\$hei,cex=1.1, col=transp("black",0.85))
points(x=roxy\$pos, y=roxy\$hei, pch=20, col=transp("red",0.7))
legend("topleft", legend=paste("Correlation found:",round(max(cors), digits=4), sep=""), bty="n")
legend("topright", legend=c("Peaks selected"), col=c("red"), bty = "n", pch=c(20), cex=0.85)
}
}else{
#### get the length of the peak region
le <- length(roxy\$pos[1]:(roxy\$pos[length(roxy\$pos)]))
#### create a vector of expected indexes according to our ladder
#### number of possible tries
tries <- length(seq(tra[length(tra)], roxy\$pos[length(roxy\$pos)], 1))+100

# define function returning absolute values
abso <- function(test1, pos){
test1 <- matrix(test1, nrow=1)
res <- apply(test1,2,function(x, y){
xx1 <- abs(as.vector(x) - y)
z <- y[which(xx1 == min(xx1))[1]]
return(z)}, y=pos)
return(res)
}

# get all possible tests
all.tests <- apply(matrix(1:tries,nrow=1), 2, function(q1, q2){q1+q2},q2=tra)
# get all possible absolute differences
all.abso <- apply(all.tests, 2, abso, pos=roxy\$pos)
# get all possible correlations
all.cors <- apply(all.abso, 2, function(x,y){cor(x,y)}, y=ladder)
# get all sum of squares
step1 <- (all.tests - all.abso)^2
all.ss2 <- apply(step1, 2, sum)
# get all standarized variances and final response
sss2 <- abs ( (all.ss2 - mean(all.ss2) )/ sd(all.ss2)) # the smaller the better
response <- all.cors/sss2

vv4 <- which(response >= min(sort(response, decreasing=TRUE)[1:5]))
#vv4 <- which(response >= .99)
reduced <- sort(unique(as.vector(all.tests[,vv4])))
reduced2 <- roxy\$pos[which(roxy\$pos >= min(reduced))]

## if there's still just toomany roxy after the reduced search
if(length(reduced2) >= nnn & length(reduced2) < nnn+8){
cors <- apply(all.combs, 2, function(x,y){cor(x,y)}, y=ladder)
# positions of roxy found
found <- all.combs[,which(cors == max(cors))]
if(draw == TRUE){
#facto <- roundUP(max(roxy\$hei) / 10)
#yyline <- seq(0,2*max(roxy\$hei), by=facto)

#if((length(x) - min(roxy\$pos)) > 500){nono <- min(roxy\$pos) - 200}else{nono <- 0}
plot(x, type="l", xaxt="n", ylim=c(0,(max(roxy\$hei)+1000)), cex.axis=0.6, las=2, xlim=c((min(roxy\$pos)-100),(max(roxy\$pos)+100)), col=transp("grey35",0.7), ylab="RFU", xlab="", lwd=2, main=attributes(x)\$mycomm, cex.main=cex.title)
#rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col = "gray90")

#abline(h=yyline, col="white", lwd=1.5, lty=3)
#lines(x=x, lwd=2, col=transp("grey25",0.6))
#polygon(x=x, col=transp("grey35",0.6))
axis(1, at=roxy\$pos, labels=roxy\$wei, cex.axis=0.6)
#abline(v=found, col="red", lty=3)
points(x=roxy\$pos, y=roxy\$hei,cex=1.1, col=transp("black",0.85))
points(x=roxy\$pos, y=roxy\$hei, pch=20, col=transp("red",0.7))
legend("topleft", legend=paste("Correlation found:",round(max(cors), digits=4), sep=""), bty="n")
legend("topright", legend=c("Peaks selected"), col=c("red"), bty = "n", pch=c(20), cex=0.85)

}
#which(sss2 == min(sss2))
# reduced search to maximum correlations
}else{# there's actually no ladder so no good correlation was found
roxy <- list(pos=seq(1,nnn) + rnorm(nnn,0,1), hei=seq(1,nnn)+ rnorm(nnn,0,1), wei= ladder)
print("Friend I don't think you have ladder in this sample or you used the wrong ladder, look the plot")
plot(x, type="l", main=attributes(x)\$mycomm, cex.main=cex.title)
}
}
}
### end of method "red"
#####################
#####################
#####################
if(method == "cor"){ #EXHAUSTIVE CORRELATION
print(paste("WOOW too many peaks in this",who,"!! low thresholds throw too many noisy peaks, consider increasing the initial thresold for your ladder, the number of possible combinations is too high to be computed, we will have to do 15,000 samples and get the most likely sizing, you better double check this sample"))

#thresh = init.thresh
#roxy <- big.peaks.col(x, thresh)
nono <- which(roxy\$pos < avoid)
roxy <- list(pos = roxy\$pos[ -nono], hei = roxy\$hei[-nono])
roxy <- separate(roxy, shift=sep.index, type="pos")
for(k in 1:15000){pos.mod[,k] <- sort(sample(roxy\$pos, size=length(ladder), replace=FALSE), decreasing=FALSE)}
v <- which(dd == max(dd))[1]
v2 <- which(roxy\$pos %in% pos.mod[, v])
roxy <- list(pos = pos.mod[, v], hei = roxy\$hei[v2], wei = ladder)
}else{
if(mi > length(roxy\$pos)){
print(paste("ERROR!! using the initial threshold you specified we did not find enough peaks for", who,", please try reducing the 'init.thresh' or check the plot for this plant using 'detect.ladder' function, maybe this sample did not have ladder :)"))
}else{
pos.mod <- combn(roxy\$pos, m=mi)
v <- which(dd == max(dd))
v2 <- which(roxy\$pos %in% pos.mod[,v])
}
}
}
# END OF METHOD "cor"
#######################
#####################
#####################
#####################
if(method == "ci"){
z <- which(roxy\$hei ==  max(roxy\$hei))
mm <- median(roxy\$hei[-z]) # get the median height of the peaks
se2.low <- (sd(roxy\$hei[-z])/sqrt(length(roxy\$hei[-z]))) * ci.low # produce the confidnce interval
se2.upp <- (sd(roxy\$hei[-z])/sqrt(length(roxy\$hei[-z]))) * ci.upp
v <- which( (roxy\$hei > (mm-se2.low)) & (roxy\$hei < (mm+se2.upp))) # reduce the ladder to the peaks inside the 95% confidence interval
roxy <- list(pos=roxy\$pos[v], hei=roxy\$hei[v])
## keep looking
vv <- which(diff(roxy\$pos) < dev)
vv2 <- vv + 1
# start condition
while(length(vv) > 0){
keep <- numeric()
for(h in 1:length(vv)){
a1 <- vv[h]
a2 <- vv2[h]
a3 <- c(roxy\$hei[a1],roxy\$hei[a2])
a4 <- c(a1,a2)
keep[h] <- (a4[which(a3 == max(a3))])[1]
}
keep <- unique(keep)
'%!in%' <- function(x,y)!('%in%'(x,y))
keep2 <- unique(c(vv,vv2)[which(c(vv,vv2) %!in% keep)])
roxy <- list(pos=roxy\$pos[-keep2], hei=roxy\$hei[-keep2])
# check again
vv <- which(diff(roxy\$pos) < dev)
vv2 <- vv + 1
}

s1 <- length(roxy\$pos)- (length(ladder) - 1)
s2 <- length(roxy\$pos)
if(s1 <= 0){
if(warn==TRUE){
print("Are you sure this is a ladder channel, we did not find a clear pattern, stop a minute to check the plot")
}
if(draw != F){
plot(x, ylim=c(0,mm+se2.upp), type="l", main=attributes(x)\$mycomm, cex.main=cex.title)
}
#roxy <- list(pos=roxy\$pos, hei=roxy\$hei, wei= ladder)
}else{roxy <- list(pos=roxy\$pos[s1:s2], hei=roxy\$hei[s1:s2], wei= ladder)}
}
################
##############################################################
if(method == "cor" | method == "ci"){
if(draw == TRUE){
xx <- roxy\$wei
yy <- roxy\$pos
mod <- lm(yy~xx)

plot(x, type="l", main=attributes(x)\$mycomm,cex.main=cex.title, xlab="", ylab="Intensity", xaxt="n", col=transp("grey39",0.6), xlim=c(0,(max(roxy\$pos)+100)), ylim=c(-100,max(roxy\$hei)))
if(method == "cor"){
legend("topleft", legend=paste("Correlation=",round(max(dd)), sep=""), bty="n")
}
abline(v=roxy[[1]], col="red", lty=3)
if(method == "ci"){
abline(h=mm, col="blue", lty=3)
abline(h=(mm+se2.upp), col="blue", lty=3)
abline(h=(mm-se2.low), col="blue", lty=3)
legend("topright", legend=c("90% CI", "Peaks selected"), col=c("blue", "red"), bty = "n", lty=c(3,3), cex=1)
}else{legend("topright", legend=c("Peaks selected"), col=c("red"), bty = "n", lty=c(3), cex=1)}