system("ms 50 10 -T -I 3 20 20 10 -n 3 4 -ma x 0.2 8 10 x 0 0.3 4 x -seed 355 339 433 | tail -n +4 | grep -v // >treefile") system("seq-gen -mHKY -l 1000 -s .2 sequences2.phy") outfile <- "simulated.mig" nloci <- 10 samplesizes <- c(20,20,10) seqnames <- array("",dim=c(nloci,sum(samplesizes))) seqstrings <- array("",dim=c(nloci,sum(samplesizes))) seq <- readLines("sequences2.phy") L <- 1 for(i in 1:nloci) { locisize <- as.numeric(strsplit(seq[L]," ")[[1]][2]) # Number of sequences for this locus (sould be sum(samplesizes)) seqs <- strsplit(seq[(L+1):(L+locisize)]," +") seqs <- t(array(unlist(seqs),c(2,locisize))) ord <- order(as.numeric(seqs[,1])) seqnames[i,]<-paste(seqs[ord,1],sapply(lapply(10-nchar(seqs[ord,1]),rep,x=" "),paste,sep="",collapse=""),sep="") seqstrings[i,]<-seqs[ord,2] L<-L+locisize+1 } cat("3 10 simulated data with 3 subpopulations and 10 loci","\n",file=outfile) cat(rep(1000,nloci),"\n",file=outfile,append=TRUE) cat(rep(samplesizes[1],nloci), "West","\n",file=outfile,append=TRUE) for(i in 1:nloci) cat(paste(seqnames[i,1:samplesizes[1]],seqstrings[i,1:samplesizes[1]]),sep="\n",file=outfile,append=TRUE) cat(rep(samplesizes[2],nloci), "South","\n",file=outfile,append=TRUE) for(i in 1:nloci) cat(paste(seqnames[i,(samplesizes[1]+1):sum(samplesizes[1:2])],seqstrings[i,(samplesizes[1]+1):sum(samplesizes[1:2])]),sep="\n",file=outfile,append=TRUE) cat(rep(samplesizes[3],nloci), "South","\n",file=outfile,append=TRUE) for(i in 1:nloci) cat(paste(seqnames[i,(sum(samplesizes[1:2])+1):sum(samplesizes)],seqstrings[i,(sum(samplesizes[1:2])+1):sum(samplesizes)]),sep="\n",file=outfile,append=TRUE)