.Rhistory

<<<<<<< HEAD
trace$G[summed==0]<-0
return(trace)
}
trace_1<-get_ratio(trace_1)
trace_2<-get_ratio(trace_2)
trace_3<-get_ratio(trace_3)
scale_optim<-function(par, trace_1, trace_2, trace_3){
#if(all(par[1]>0, par[1]<1, par[2]>0, par[2]<1)){
scale_1<-par[1]
scale_2<-par[2]
ssq<-sum(((scale_1*trace_1+scale_2*trace_2)-trace_3)^2)
return(ssq)#}else{
#return(NA)
# }
}
scale_fit<-function(trace_1,trace_2,trace_3){
# #Extract base trace information
trace_1<-trace_1[,names(trace_1) %in% c("A","T", "C","G")]
trace_2<-trace_2[,names(trace_2) %in% c("A","T", "C","G")]
trace_3<-trace_3[,names(trace_3) %in% c("A","T", "C","G")]
=======
trace_2=trace_2,
trace_3=trace_3, lower=c(0,0), upper=c(1,1)))
?DEopt
system.time(scale.fit<-DEopt(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, lower=c(0,0), upper=c(1,1)))
OF <- scale_optim               ### see ?testFunctions
algo <- list(nP = 50L,          ### population size
nG = 300L,         ### number of generations
F = 0.01,          ### step size
CR = 0.9,          ### prob of crossover
min = c(0, 0),  ### range for initial population
max = c(1, 1))
system.time(scale.fit<-DEopt(OF=scale_optim,
algo,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3))
OF <- scale_optim               ### see ?testFunctions
algo <- list(nP = 10L,          ### population size
nG = 25L,         ### number of generations
F = 0.01,          ### step size
CR = 0.9,          ### prob of crossover
min = c(0, 0),  ### range for initial population
max = c(1, 1))
system.time(scale.fit<-DEopt(OF=scale_optim,
algo,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3))
?optim
?PSopt
system.time(scale.fit<-PSopt(OF=scale_optim,
algo,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3))
algo <- list(nP = 50L,          ### population size
nG = 300L,         ### number of generations
F = 0.01,          ### step size
CR = 0.9,          ### prob of crossover
min = c(0, 0),  ### range for initial population
max = c(1, 1))
system.time(scale.fit<-DEopt(OF=scale_optim,
algo,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3))
algo <- list(nP = 50L,          ### population size
nG = 25L,         ### number of generations
F = 0.01,          ### step size
CR = 0.9,          ### prob of crossover
min = c(0, 0),  ### range for initial population
max = c(1, 1))
system.time(scale.fit<-DEopt(OF=scale_optim,
algo,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3))
algo <- list(nP = 25L,          ### population size
nG = 25L,         ### number of generations
F = 0.01,          ### step size
CR = 0.9,          ### prob of crossover
min = c(0, 0),  ### range for initial population
max = c(1, 1))
system.time(scale.fit<-DEopt(OF=scale_optim,
algo,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3))
system.time(scale.fit<-PSopt(OF=scale_optim,
algo,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3))
?optim
?optim
system.time(scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=500))
)
scale_optim<-function(par, trace_1, trace_2, trace_3){
scale_1<-par[1]
scale_2<-par[2]
ssq<-sum(((scale_1*trace_1+scale_2*trace_2)-trace_3)^2)
return(ssq)
}
system.time(scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=500)))
system.time(scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=50)))
scale.fit
system.time(scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=10)))
length(-150:150)
length(-150:150)/31
31^2/11^2
scale.fit
optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=100))
optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=100))
scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=5))
scale.fit
optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=5, ndeps=c(0.1,0.1)))
optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=10, ndeps=c(0.1,0.1)))
optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=100, ndeps=c(0.1,0.1)))
optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=20, ndeps=c(0.1,0.1)))
scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=15, ndeps=c(0.1,0.1)))
scale.fit
scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=10, ndeps=c(0.1,0.1)))
scale.fit
system.time(scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=15, ndeps=c(0.1,0.1))))
scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=15))
scale.fit
scale_optim<-function(par, trace_1, trace_2, trace_3){
if(all(par[1]>0, par[1]<1, par[2]>0, par[2]<1)){
scale_1<-par[1]
scale_2<-par[2]
ssq<-sum(((scale_1*trace_1+scale_2*trace_2)-trace_3)^2)
return(ssq)}else{
return(NA)
}
}
system.time(scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3, control=list(maxit=15)))
system.time(scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3,
method="L-BFGS-B",
lower=c(0,0),
upper=c(1,1)))
scale_optim<-function(par, trace_1, trace_2, trace_3){
#if(all(par[1]>0, par[1]<1, par[2]>0, par[2]<1)){
scale_1<-par[1]
scale_2<-par[2]
ssq<-sum(((scale_1*trace_1+scale_2*trace_2)-trace_3)^2)
return(ssq)#}else{
#return(NA)
# }
}
system.time(scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3,
method="L-BFGS-B",
lower=c(0,0),
upper=c(1,1)))
scale.fit
>>>>>>> no change
scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3,
method="L-BFGS-B",
lower=c(0,0),
upper=c(1,1),
<<<<<<< HEAD
control=list(maxit=50))
return(scale.fit)}
scale.fit<-scale_fit(trace_1,trace_2,trace_3)
fit.com<-scale.fit$par[1]*trace_1[,names(trace_1) %in% c("A","T", "C","G")]+scale.fit$par[2]*trace_2[,names(trace_1) %in% c("A","T", "C","G")]
trace_1<-temp_1
trace_2<-temp_2
trace_3<-temp_3
plot_trace(fit.com,
region=c(5000,10000),
print.plot=F,
main="trace_3",
split=F,
base.call=F,
ylim=c(0,2050))
fit.com$datum<-1:nrow(fit.com)
plot_trace(fit.com,
region=c(5000,10000),
print.plot=F,
main="trace_3",
split=F,
base.call=F,
ylim=c(0,2050))
fit.com$datum<-1:nrow(fit.com)
plot_trace(fit.com,
#region=c(5000,10000),
print.plot=F,
main="trace_3",
split=F,
base.call=F,
ylim=c(0,2050))
head(fit.com)
trace_1<-temp_1
trace_2<-temp_2
trace_3<-temp_3
scale.fit<-scale_fit(trace_1,trace_2,trace_3)
fit.com<-scale.fit$par[1]*trace_1[,names(trace_1) %in% c("A","T", "C","G")]+scale.fit$par[2]*trace_2[,names(trace_1) %in% c("A","T", "C","G")]
fit.com$datum<-1:nrow(fit.com)
plot_trace(fit.com,
#region=c(5000,10000),
print.plot=F,
main="trace_3",
split=F,
base.call=F,
ylim=c(0,2050))
trace_1<-adam_test_data[[1]]
trace_2<-adam_test_data[[2]]
trace_3<-adam_test_data[[3]]
#Scale two visualy aligned traces
#Save a temp of traces
temp_1<-trace_1
temp_2<-trace_2
temp_3<-trace_3
#Get visually aligned region
# trace_1<-trace_1[5000:5500,]
# trace_2<-trace_2[4955:5455,]
# trace_3<-trace_3[4910:5410,]
# trace_1<-trace_1[6000:6500,]
# trace_2<-trace_2[5955:6455,]
# trace_3<-trace_3[5890:6390,]
trace_1<-trace_1[6000:6150,]
trace_2<-trace_2[5955:6105,]
trace_3<-trace_3[5890:6040,]
scale_fit<-function(trace_1,trace_2,trace_3){
#Remove minimum values
rm.low<-colwise(function(x){
return(ifelse(x < 500, 0, x))})
rm.low.df<-function(trace){
trace[,names(trace) %in% c("A","T","C","G")]<-rm.low(trace[,names(trace) %in% c("A","T","C","G")])
return(trace)
}
trace_1<-rm.low.df(trace_1)
trace_2<-rm.low.df(trace_2)
trace_3<-rm.low.df(trace_3)
####
# Get ratios
###
get_ratio<-function(trace){
summed<-apply(trace[,names(trace) %in% c("A","T","C","G")], MARGIN=1, FUN=sum)
trace[,names(trace) %in% c("A","T","C","G")]<-trace[,names(trace) %in% c("A","T","C","G")]/summed
trace$A[summed==0]<-0
trace$T[summed==0]<-0
trace$C[summed==0]<-0
trace$G[summed==0]<-0
return(trace)
}
trace_1<-get_ratio(trace_1)
trace_2<-get_ratio(trace_2)
trace_3<-get_ratio(trace_3)
scale_optim<-function(par, trace_1, trace_2, trace_3){
#if(all(par[1]>0, par[1]<1, par[2]>0, par[2]<1)){
scale_1<-par[1]
scale_2<-par[2]
ssq<-sum(((scale_1*trace_1+scale_2*trace_2)-trace_3)^2)
return(ssq)#}else{
#return(NA)
# }
}
# #Extract base trace information
trace_1<-trace_1[,names(trace_1) %in% c("A","T", "C","G")]
trace_2<-trace_2[,names(trace_2) %in% c("A","T", "C","G")]
trace_3<-trace_3[,names(trace_3) %in% c("A","T", "C","G")]
=======
control=list(maxit=15))
#Align and scale two traces so they add up to a third
#Align using all base traces (A and B and C...)
source("./R/pad_shift.R")
source("./R/blank_length.R")
align_trace<-function(trace_1, trace_2, trace_3, trace_values=F){
if(!require(NMOF)){install.packages("NMOF")}
if(!require(multicore)){install.packages("multicore")}
if(!require(parallel)){install.packages("parallel")}
temp_1<-trace_1
temp_2<-trace_2
temp_3<-trace_3
#Extract base trace information
trace_1<-trace_1[,names(trace_1) %in% c("A","T", "C","G")]
trace_2<-trace_2[,names(trace_2) %in% c("A","T", "C","G")]
trace_3<-trace_3[,names(trace_3) %in% c("A","T", "C","G")]
#Pad so that all traces have the same length
lengthened<-blank_length(trace_1, trace_2, trace_3)
trace_1<-lengthened[[1]]
trace_2<-lengthened[[2]]
trace_3<-lengthened[[3]]
#Define a scale optimization function
scale_optim<-function(par, trace_1, trace_2, trace_3){
#if(all(par[1]>0, par[1]<1, par[2]>0, par[2]<1)){
scale_1<-par[1]
scale_2<-par[2]
ssq<-sum(((scale_1*trace_1+scale_2*trace_2)-trace_3)^2)
return(ssq)#}else{
#return(NA)
# }
}
#############################
#Define a shift optimziation function that calls the scale optim internaly
shift_optim<-function(par, trace_1, trace_2, trace_3){
shift_1<-par[1]
shift_2<-par[2]
#############################
#Pad with assigned shifts
#Shift trace 1 as function of 2 and 3
padded_1<-pad.shift(shift_1, trace_1, trace_2, trace_3)
trace_1<-padded_1[[1]]
trace_2<-padded_1[[2]]
trace_3<-padded_1[[3]]
#Shift trace 2 as a function of 1 and 3
padded_2<-pad.shift(shift_2, trace_2, trace_1, trace_3)
trace_1<-padded_2[[1]]
trace_2<-padded_2[[2]]
trace_3<-padded_2[[3]]
#Fit the scale to the optimal alignement
scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3,
method="L-BFGS-B",
lower=c(0,0),
upper=c(1,1),
control=list(maxit=15))
return(scale.fit$value[1])
}
fit<-gridSearch(fun=shift_optim,
lower=-150,
upper=150,
npar=2,
n=11,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3,
method="multicore",
mc.control = list(mc.cores=3))
shifts<-fit$minlevels
shift_1<-shifts[1]
shift_2<-shifts[2]
#Use optimal values to create new trace vectore
#Shift trace 1 as function of 2 and 3
padded_1<-pad.shift(shift_1, trace_1, trace_2, trace_3)
trace_1<-padded_1[[1]]
trace_2<-padded_1[[2]]
trace_3<-padded_1[[3]]
#Shift trace 2 as a function of 1 and 3
padded_2<-pad.shift(shift_2, trace_2, trace_1, trace_3)
trace_1<-padded_2[[1]]
trace_2<-padded_2[[2]]
trace_3<-padded_2[[3]]
#Obtain the scale parameter on the best aligned sequences
scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3)
scale_1<-scale.fit$par[1]
scale_2<-scale.fit$par[2]
#Do it using an Differential evolution algorithm
# algo <- list(nP = 50L,          ### population size
#              nG = 300L,         ### number of generations
#              F = 0.01,          ### step size
#              CR = 0.9,          ### prob of crossover
#              min = c(0, 0),  ### range for initial population
#              max = c(1, 1))
#
# system.time(scale.fit<-DEopt(OF=scale_optim,
#                              algo,
#                                  trace_1=trace_1,
#                                  trace_2=trace_2,
#                                  trace_3=trace_3))
#
# system.time(scale.fit<-PSopt(OF=scale_optim,
#                              algo,
#                                  trace_1=trace_1,
#                                  trace_2=trace_2,
#                                  trace_3=trace_3))
shift_n_scale<-data.frame(shift_1=shift_1,
shift_2=shift_2,
scale_1=scale_1,
scale_2=scale_2)
if(trace_values){
#Scale
trace_1<-scale_1*trace_1
trace_2<-scale_2*trace_2
#Add initial data
trace_1<-join(trace_1, temp_1,
type="full",
sort=F)
trace_2<-merge(trace_2, temp_2,
all=T,
sort=F)
trace_3<-merge(trace_3, temp_3,
all=T,
sort=F)
traces<-list(trace_1,
trace_2,
trace_3)
print(shift_n_scale)
return(traces)
}else{
return(shift_n_scale)
}
}
#Test
system.time(shift<-align_trace(trace_1, trace_2, trace_3))
#Optimization effort scales with the square (two parameter) of the number of steps
#
shift
41^2/11^2*121
41^2/11^2*121/60
81^2/11^2*121/60
81^2/11^2*121/60/60
400/80
#Align and scale two traces so they add up to a third
#Align using all base traces (A and B and C...)
source("./R/pad_shift.R")
source("./R/blank_length.R")
align_trace<-function(trace_1, trace_2, trace_3, trace_values=F){
if(!require(NMOF)){install.packages("NMOF")}
if(!require(multicore)){install.packages("multicore")}
if(!require(parallel)){install.packages("parallel")}
temp_1<-trace_1
temp_2<-trace_2
temp_3<-trace_3
#Extract base trace information
trace_1<-trace_1[,names(trace_1) %in% c("A","T", "C","G")]
trace_2<-trace_2[,names(trace_2) %in% c("A","T", "C","G")]
trace_3<-trace_3[,names(trace_3) %in% c("A","T", "C","G")]
#Pad so that all traces have the same length
lengthened<-blank_length(trace_1, trace_2, trace_3)
trace_1<-lengthened[[1]]
trace_2<-lengthened[[2]]
trace_3<-lengthened[[3]]
#Define a scale optimization function
scale_optim<-function(par, trace_1, trace_2, trace_3){
#if(all(par[1]>0, par[1]<1, par[2]>0, par[2]<1)){
scale_1<-par[1]
scale_2<-par[2]
ssq<-sum(((scale_1*trace_1+scale_2*trace_2)-trace_3)^2)
return(ssq)#}else{
#return(NA)
# }
}
#############################
#Define a shift optimziation function that calls the scale optim internaly
shift_optim<-function(par, trace_1, trace_2, trace_3){
shift_1<-par[1]
shift_2<-par[2]
#############################
#Pad with assigned shifts
#Shift trace 1 as function of 2 and 3
padded_1<-pad.shift(shift_1, trace_1, trace_2, trace_3)
trace_1<-padded_1[[1]]
trace_2<-padded_1[[2]]
trace_3<-padded_1[[3]]
#Shift trace 2 as a function of 1 and 3
padded_2<-pad.shift(shift_2, trace_2, trace_1, trace_3)
trace_1<-padded_2[[1]]
trace_2<-padded_2[[2]]
trace_3<-padded_2[[3]]
#Fit the scale to the optimal alignement
>>>>>>> no change
scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3,
method="L-BFGS-B",
lower=c(0,0),
upper=c(1,1),
<<<<<<< HEAD
control=list(maxit=50))
return(scale.fit)}
scale.fit<-scale_fit(trace_1,trace_2,trace_3)
fit.com<-scale.fit$par[1]*trace_1[,names(trace_1) %in% c("A","T", "C","G")]+scale.fit$par[2]*trace_2[,names(trace_1) %in% c("A","T", "C","G")]
fit.com$datum<-1:nrow(fit.com)
plot_trace(fit.com,
#region=c(5000,10000),
print.plot=F,
main="trace_3",
split=F,
base.call=F,
ylim=c(0,2050))
trace_1<-temp_1
trace_2<-temp_2
trace_3<-temp_3
#Scale two visualy aligned traces
#Save a temp of traces
temp_1<-trace_1
temp_2<-trace_2
temp_3<-trace_3
#Get visually aligned region
# trace_1<-trace_1[5000:5500,]
# trace_2<-trace_2[4955:5455,]
# trace_3<-trace_3[4910:5410,]
# trace_1<-trace_1[6000:6500,]
# trace_2<-trace_2[5955:6455,]
# trace_3<-trace_3[5890:6390,]
trace_1<-trace_1[6000:6150,]
trace_2<-trace_2[5955:6105,]
trace_3<-trace_3[5890:6040,]
scale_fit<-function(trace_1,trace_2,trace_3){
#Remove minimum values
rm.low<-colwise(function(x){
return(ifelse(x < 500, 0, x))})
rm.low.df<-function(trace){
trace[,names(trace) %in% c("A","T","C","G")]<-rm.low(trace[,names(trace) %in% c("A","T","C","G")])
return(trace)
}
trace_1<-rm.low.df(trace_1)
trace_2<-rm.low.df(trace_2)
trace_3<-rm.low.df(trace_3)
####
# Get ratios
###
get_ratio<-function(trace){
summed<-apply(trace[,names(trace) %in% c("A","T","C","G")], MARGIN=1, FUN=sum)
trace[,names(trace) %in% c("A","T","C","G")]<-trace[,names(trace) %in% c("A","T","C","G")]/summed
trace$A[summed==0]<-0
trace$T[summed==0]<-0
trace$C[summed==0]<-0
trace$G[summed==0]<-0
return(trace)
}
trace_1<-get_ratio(trace_1)
trace_2<-get_ratio(trace_2)
trace_3<-get_ratio(trace_3)
scale_optim<-function(par, trace_1, trace_2, trace_3){
#if(all(par[1]>0, par[1]<1, par[2]>0, par[2]<1)){
scale_1<-par[1]
scale_2<-par[2]
ssq<-sum(((scale_1*trace_1+scale_2*trace_2)-trace_3)^2)
return(ssq)#}else{
#return(NA)
# }
}
# #Extract base trace information
trace_1<-trace_1[,names(trace_1) %in% c("A","T", "C","G")]
trace_2<-trace_2[,names(trace_2) %in% c("A","T", "C","G")]
trace_3<-trace_3[,names(trace_3) %in% c("A","T", "C","G")]
scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3,
method="L-BFGS-B",
lower=c(0,0),
upper=c(1,1),
control=list(maxit=50))
return(scale.fit)}
scale.fit<-scale_fit(trace_1,trace_2,trace_3)
fit.com<-scale.fit$par[1]*trace_1[,names(trace_1) %in% c("A","T", "C","G")]+scale.fit$par[2]*trace_2[,names(trace_1) %in% c("A","T", "C","G")]
fit.com$datum<-1:nrow(fit.com)
plot_trace(fit.com,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F,
ylim=c(0,2050))
trace_1<-temp_1
trace_2<-temp_2
trace_3<-temp_3
trace_1<-trace_1[6000:6150,]
trace_2<-trace_2[5955:6105,]
trace_3<-trace_3[5890:6040,]
#Remove minimum values
rm.low<-colwise(function(x){
return(ifelse(x < 500, 0, x))})
rm.low.df<-function(trace){
trace[,names(trace) %in% c("A","T","C","G")]<-rm.low(trace[,names(trace) %in% c("A","T","C","G")])
return(trace)
}
trace_1<-rm.low.df(trace_1)
trace_2<-rm.low.df(trace_2)
trace_3<-rm.low.df(trace_3)
####
# Get ratios
###
get_ratio<-function(trace){
summed<-apply(trace[,names(trace) %in% c("A","T","C","G")], MARGIN=1, FUN=sum)
trace[,names(trace) %in% c("A","T","C","G")]<-trace[,names(trace) %in% c("A","T","C","G")]/summed
trace$A[summed==0]<-0
trace$T[summed==0]<-0
trace$C[summed==0]<-0
trace$G[summed==0]<-0
return(trace)
}
trace_1<-get_ratio(trace_1)
trace_2<-get_ratio(trace_2)
trace_3<-get_ratio(trace_3)
scale_optim<-function(par, trace_1, trace_2, trace_3){
#if(all(par[1]>0, par[1]<1, par[2]>0, par[2]<1)){
scale_1<-par[1]
scale_2<-par[2]
ssq<-sum(((scale_1*trace_1+scale_2*trace_2)-trace_3)^2)
return(ssq)#}else{
#return(NA)
# }
}
# #Extract base trace information
trace_1<-trace_1[,names(trace_1) %in% c("A","T", "C","G")]
trace_2<-trace_2[,names(trace_2) %in% c("A","T", "C","G")]
trace_3<-trace_3[,names(trace_3) %in% c("A","T", "C","G")]
trace_1<-temp_1
trace_2<-temp_2
trace_3<-temp_3
#Remove minimum values
rm.low<-colwise(function(x){
return(ifelse(x < 500, 0, x))})
rm.low.df<-function(trace){
trace[,names(trace) %in% c("A","T","C","G")]<-rm.low(trace[,names(trace) %in% c("A","T","C","G")])
return(trace)
}
trace_1<-rm.low.df(trace_1)
trace_2<-rm.low.df(trace_2)
trace_3<-rm.low.df(trace_3)
####
# Get ratios
###
get_ratio<-function(trace){
summed<-apply(trace[,names(trace) %in% c("A","T","C","G")], MARGIN=1, FUN=sum)
trace[,names(trace) %in% c("A","T","C","G")]<-trace[,names(trace) %in% c("A","T","C","G")]/summed
trace$A[summed==0]<-0
trace$T[summed==0]<-0
trace$C[summed==0]<-0
trace$G[summed==0]<-0
return(trace)
}
trace_1<-get_ratio(trace_1)
trace_2<-get_ratio(trace_2)
trace_3<-get_ratio(trace_3)
scale_optim<-function(par, trace_1, trace_2, trace_3){
#if(all(par[1]>0, par[1]<1, par[2]>0, par[2]<1)){
scale_1<-par[1]
scale_2<-par[2]
ssq<-sum(((scale_1*trace_1+scale_2*trace_2)-trace_3)^2)
return(ssq)#}else{
#return(NA)
# }
}
# #Extract base trace information
trace_1<-trace_1[,names(trace_1) %in% c("A","T", "C","G")]
trace_2<-trace_2[,names(trace_2) %in% c("A","T", "C","G")]
trace_3<-trace_3[,names(trace_3) %in% c("A","T", "C","G")]
fit.com<-scale.fit$par[1]*trace_1[,names(trace_1) %in% c("A","T", "C","G")]+scale.fit$par[2]*trace_2[,names(trace_1) %in% c("A","T", "C","G")]
fit.com$datum<-1:nrow(fit.com)
plot_trace(fit.com,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F,
ylim=c(0,2050))
plot_trace(fit.com,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F)
trace_1<-adam_test_data[[1]]
trace_2<-adam_test_data[[2]]
trace_3<-adam_test_data[[3]]
trace_1<-trace_1[6000:6150,]
trace_2<-trace_2[5955:6105,]
trace_3<-trace_3[5890:6040,]
#Remove minimum values
rm.low<-colwise(function(x){
return(ifelse(x < 500, 0, x))})
rm.low.df<-function(trace){
trace[,names(trace) %in% c("A","T","C","G")]<-rm.low(trace[,names(trace) %in% c("A","T","C","G")])
return(trace)
}
trace_1<-rm.low.df(trace_1)
trace_2<-rm.low.df(trace_2)
trace_3<-rm.low.df(trace_3)
####
# Get ratios
###
get_ratio<-function(trace){
summed<-apply(trace[,names(trace) %in% c("A","T","C","G")], MARGIN=1, FUN=sum)
trace[,names(trace) %in% c("A","T","C","G")]<-trace[,names(trace) %in% c("A","T","C","G")]/summed
trace$A[summed==0]<-0
trace$T[summed==0]<-0
trace$C[summed==0]<-0
trace$G[summed==0]<-0
return(trace)
}
trace_1<-get_ratio(trace_1)
trace_2<-get_ratio(trace_2)
trace_3<-get_ratio(trace_3)
scale_optim<-function(par, trace_1, trace_2, trace_3){
#if(all(par[1]>0, par[1]<1, par[2]>0, par[2]<1)){
scale_1<-par[1]
scale_2<-par[2]
ssq<-sum(((scale_1*trace_1+scale_2*trace_2)-trace_3)^2)
return(ssq)#}else{
#return(NA)
# }
}
# #Extract base trace information
trace_1<-trace_1[,names(trace_1) %in% c("A","T", "C","G")]
trace_2<-trace_2[,names(trace_2) %in% c("A","T", "C","G")]
trace_3<-trace_3[,names(trace_3) %in% c("A","T", "C","G")]
scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3,
method="L-BFGS-B",
lower=c(0,0),
upper=c(1,1),
control=list(maxit=50))
scale.fit<-scale_fit(trace_1,trace_2,trace_3)
fit.com<-scale.fit$par[1]*trace_1[,names(trace_1) %in% c("A","T", "C","G")]+scale.fit$par[2]*trace_2[,names(trace_1) %in% c("A","T", "C","G")]
fit.com$datum<-1:nrow(fit.com)
plot_trace(fit.com,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F)
trace_1$datum<-1:nrow(trace_1)
trace_2$datum<-1:nrow(trace_2)
trace_3$datum<-1:nrow(trace_3)
plot_trace(trace_1,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F)
plot_trace(trace_1,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F)+opts(legend.position="none")
plot_trace(trace_2,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F)+opts(legend.position="none")
plot_trace(fit.com,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F)
plot_trace(trace_3,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F)+opts(legend.position="none")
plot_trace(trace_1,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F)+opts(legend.position="none")
plot_trace(trace_1,
#region=c(5000,10000),
print.plot=F,
main="trace 1",
split=F,
base.call=F)+opts(legend.position="none")
plot_trace(trace_2,
#region=c(5000,10000),
print.plot=F,
main="trace 2",
split=F,
base.call=F)+opts(legend.position="none")
plot_trace(fit.com,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F)
plot_trace(trace_3,
#region=c(5000,10000),
print.plot=F,
main="Community trace",
split=F,
base.call=F)+opts(legend.position="none")
plot_trace(fit.com,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F)+opts(legend.position="none")
plot_trace(trace_3,
#region=c(5000,10000),
print.plot=F,
main="Community trace",
split=F,
base.call=F)+opts(legend.position="none")
trace_1<-temp_1
trace_2<-temp_2
trace_3<-temp_3
trace_1<-trace_1[6000:6150,]
trace_2<-trace_2[5955:6105,]
trace_3<-trace_3[5890:6040,]
fit.com<-scale.fit$par[1]*trace_1[,names(trace_1) %in% c("A","T", "C","G")]+scale.fit$par[2]*trace_2[,names(trace_1) %in% c("A","T", "C","G")]
plot_trace(trace_3,
#region=c(5000,10000),
print.plot=F,
main="Community trace",
split=F,
base.call=F)+opts(legend.position="none")
plot_trace(fit.com,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F)+opts(legend.position="none")
fit.com$datum<-1:nrow(fit.com)
plot_trace(fit.com,
#region=c(5000,10000),
print.plot=F,
main="Fit trace",
split=F,
base.call=F)+opts(legend.position="none")
test<-data.frame(x=1:12, y=c(0,1,0,1,0,0,1,0,0,0,1,0))
qplot(data=test, x=x, y=y, geom="smooth")
test<-data.frame(y=c(0,1,0,1,0,0,1,0,0,0,1,0))
zero.pad<-5
spaced.trace<-data.frame(y=0)
for(i in 1:nrow(test)){
temp<-rbind(test[i,],
data.frame(y=rep(0,zero.pad)))
spaced.trace<-rbind(spaced.trace, temp)}
spaced.trace
spaced.trace<-smooth_trace(spaced.trace, type="Savitsky-Golay")
qplot(data=spaced.trace, x=1:nrow(spaced.trace), y=y)
qplot(data=data.frame(spaced.trace), x=1:nrow(spaced.trace), y=y)
head(spaced.trace)
spaced.trace<-data.frame(spaced.trace)
qplot(data=spaced.trace, x=1:nrow(spaced.trace), y=y)
nrow(spaced.trace)
spaced.trace$datum<-1:nrow(spaced.trace)
qplot(data=spaced.trace, x=datum, y=y)
head(spaced.trace)
qplot(data=spaced.trace, x=datum, y=spaced.trace, ylab="value")
qplot(data=spaced.trace, x=datum, y=spaced.trace, ylab="value", geom="smooth")
qplot(data=spaced.trace, x=datum, y=spaced.trace, ylab="value", geom="line")
=======
control=list(maxit=15))
return(scale.fit$value[1])
}
fit<-gridSearch(fun=shift_optim,
lower=-200,
upper=200,
npar=2,
n=81,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3,
method="multicore",
mc.control = list(mc.cores=3))
shifts<-fit$minlevels
shift_1<-shifts[1]
shift_2<-shifts[2]
#Use optimal values to create new trace vectore
#Shift trace 1 as function of 2 and 3
padded_1<-pad.shift(shift_1, trace_1, trace_2, trace_3)
trace_1<-padded_1[[1]]
trace_2<-padded_1[[2]]
trace_3<-padded_1[[3]]
#Shift trace 2 as a function of 1 and 3
padded_2<-pad.shift(shift_2, trace_2, trace_1, trace_3)
trace_1<-padded_2[[1]]
trace_2<-padded_2[[2]]
trace_3<-padded_2[[3]]
#Obtain the scale parameter on the best aligned sequences
scale.fit<-optim(par=c(0.5, 0.5),
fn=scale_optim,
trace_1=trace_1,
trace_2=trace_2,
trace_3=trace_3)
scale_1<-scale.fit$par[1]
scale_2<-scale.fit$par[2]
#Do it using an Differential evolution algorithm
# algo <- list(nP = 50L,          ### population size
#              nG = 300L,         ### number of generations
#              F = 0.01,          ### step size
#              CR = 0.9,          ### prob of crossover
#              min = c(0, 0),  ### range for initial population
#              max = c(1, 1))
#
# system.time(scale.fit<-DEopt(OF=scale_optim,
#                              algo,
#                                  trace_1=trace_1,
#                                  trace_2=trace_2,
#                                  trace_3=trace_3))
#
# system.time(scale.fit<-PSopt(OF=scale_optim,
#                              algo,
#                                  trace_1=trace_1,
#                                  trace_2=trace_2,
#                                  trace_3=trace_3))
shift_n_scale<-data.frame(shift_1=shift_1,
shift_2=shift_2,
scale_1=scale_1,
scale_2=scale_2)
if(trace_values){
#Scale
trace_1<-scale_1*trace_1
trace_2<-scale_2*trace_2
#Add initial data
trace_1<-join(trace_1, temp_1,
type="full",
sort=F)
trace_2<-merge(trace_2, temp_2,
all=T,
sort=F)
trace_3<-merge(trace_3, temp_3,
all=T,
sort=F)
traces<-list(trace_1,
trace_2,
trace_3)
print(shift_n_scale)
return(traces)
}else{
return(shift_n_scale)
}
}
#Test
system.time(shift<-align_trace(trace_1, trace_2, trace_3))
#Optimization effort scales with the square (two parameter) of the number of steps
#
shift
save(shift, file="./Outputs/Shift and scale at 5 res.RData")
6142/60/60
401^2/81^2
401^2/81^2*1.706
>>>>>>> no change