R/gn3sim.R

In DNAseqtest: Generating and Testing DNA Sequences

gn3sim <-
function(theta,seqLength,merge2){
genseq2 <-function(t1,t2,seq,Sx2,Sy2,Pix,Piy){
sequ<-matrix(0,ncol=2,nrow=length(seq))
p1<-Pt(Sx2,Pix,t1)
p2<-Pt(Sy2,Piy,t2)
for(i in 1:length(seq)){
sequ[i,1]<-sample(1:4,1,prob=p1[seq[i],])
sequ[i,2]<-sample(1:4,1,prob=p2[seq[i],])
}
return(sequ)
}
pix1<-theta[1:3]
piy1<-theta[4:6]
f00<-theta[7:9]
ss<-theta[10:15]
rho<-theta[16]
pix2<-c(pix1,1-sum(pix1))
Pix2<-diag(pix2)
piy2<-c(piy1,1-sum(piy1))
Piy2<-diag(piy2)
f0<-c(f00,1-sum(f00))
Sx1<-Smatrix(ss,pix2)
Sy1<-rho*Smatrix(ss,piy2)
height<-c(theta[17:length(theta)],1)
merge<-merge2
k<-dim(merge)[1]
h<-c(0,height)
me<-merge[k,]
m<-me
t0<-h[k+1]
cc<-0
if(merge[,1][k]>0){
cc<-merge[,1][k]
v<-matrix(0,k,1)
v[1]<-merge[k,1]
cur<-0
max<-1
while(cur<max){
cur<-cur+1
for(i in c(1,2)){
if(merge[(v[cur]),i]>0){
max<-max+1
v[max]<-merge[(v[cur]),i]
}
}
cc<-matrix(0,max,1)
for(i in 1:max){
cc[i]<-v[i]
}
}
}
seqr<-matrix(0,seqLength,k-1)
fs<-matrix(0,seqLength,k+1)
seq1<-cbind(sample(1:4,seqLength,prob=f0,replace=T))
for(i in 1:k){
iii<-k-i+1
if(iii==k){
xx<-0
if(any(merge[k,1]==cc)){
pix<-pix2
Pix<-Pix2
piy<-piy2
Piy<-Piy2
Sx2<-Sx1
Sy2<-Sy1
}
else{
piy<-pix2
Piy<-Pix2
pix<-piy2
Pix<-Piy2
Sy2<-Sx1
Sx2<-Sy1
}
}
else{
if(any(iii==cc)){
pix<-pix2
Pix<-Pix2
piy<-pix2
Piy<-Pix2
Sx2<-Sx1
Sy2<-Sx1
}
else{
piy<-piy2
Piy<-Piy2
pix<-piy2
Pix<-Piy2
Sx2<-Sy1
Sy2<-Sy1
}
}
m<-rev(sort(me))
mm<-merge[m[1],]
if(me[1]>0){
t1<-t0-h[me[1]+1]
}
else{
t1<-t0
}
if(me[2]>0){
t2<-t0-h[me[2]+1]
}
else{
t2<-t0
}
if(t1<0||t1>1){
return(0)
}
if(t2<0||t2>1){
return(0)
}
hp<-rev(order(me))
if(i!=1){
seq1<-seqr[,(k-i+1)]
}
seq2<-genseq2(t1,t2,seq1,Sx2,Sy2,Pix,Piy)
if(me[1]<0){
fs[,(abs(me[1]))]<-seq2[,1]
}
else{
seqr[,(abs(me[1]))]<-seq2[,1]
}
if(me[2]<0){
fs[,(abs(me[2]))]<-seq2[,2]
}
else{
seqr[,(abs(me[2]))]<-seq2[,2]
}
if(m[1]>0)
t0<-h[m[1]+1]
else{
return(fs)
}
me<-c(mm,m[-1])
}
}