###################################################
### chunk number 1: 
###################################################
#line 55 "vignettes/clstutils/inst/doc/refSet.Rnw"
  figdir <- 'figs'
  dir.create(figdir, showWarnings=FALSE)


###################################################
### chunk number 2: loadLibs
###################################################
#line 75 "vignettes/clstutils/inst/doc/refSet.Rnw"
library(ape)
library(lattice)
library(clst)
library(clstutils)


###################################################
### chunk number 3: 
###################################################
#line 87 "vignettes/clstutils/inst/doc/refSet.Rnw"
data(seqs)
data(seqdat)


###################################################
### chunk number 4: 
###################################################
#line 97 "vignettes/clstutils/inst/doc/refSet.Rnw"
seqdat$i <- 1:nrow(seqdat)
taxa <- split(seqdat, seqdat$tax_name)
nseqs <- sapply(taxa, nrow)
nseqs


###################################################
### chunk number 5: 
###################################################
#line 112 "vignettes/clstutils/inst/doc/refSet.Rnw"
Efaecium <- taxa[['Enterococcus faecium']]$i


###################################################
### chunk number 6: 
###################################################
#line 118 "vignettes/clstutils/inst/doc/refSet.Rnw"
dmat <- ape::dist.dna(seqs[Efaecium,], pairwise.deletion=TRUE, as.matrix=TRUE, model='raw')
summary(dmat[lower.tri(dmat)])


###################################################
### chunk number 7: 
###################################################
#line 127 "vignettes/clstutils/inst/doc/refSet.Rnw"
outliers <- clstutils::findOutliers(dmat, cutoff=0.015)
table(outliers)


###################################################
### chunk number 8: 
###################################################
#line 139 "vignettes/clstutils/inst/doc/refSet.Rnw"
with(seqdat[Efaecium,], {
  prettyTree(nj(dmat), groups=ifelse(outliers,'outlier','non-outlier'),
             X=outliers, labels=ifelse(isType,gettextf('type strain (%s)', accession),NA))
})


###################################################
### chunk number 9: 
###################################################
#line 150 "vignettes/clstutils/inst/doc/refSet.Rnw"
dmats <- lapply(taxa, function(taxon) {
  ape::dist.dna(seqs[taxon$i,], pairwise.deletion=TRUE, as.matrix=TRUE, model='raw')
})


###################################################
### chunk number 10: 
###################################################
#line 161 "vignettes/clstutils/inst/doc/refSet.Rnw"
outliers <- sapply(dmats, findOutliers, cutoff=0.015)


###################################################
### chunk number 11: 
###################################################
#line 167 "vignettes/clstutils/inst/doc/refSet.Rnw"
seqdat$outlier <- FALSE
for(x in outliers){
  seqdat[match(names(x),seqdat$seqname),'outlier'] <- x
}

with(seqdat, table(tax_name, outlier))


###################################################
### chunk number 12: 
###################################################
#line 182 "vignettes/clstutils/inst/doc/refSet.Rnw"
lowerTriangle <- function(mat){mat[lower.tri(mat)]}

dists <- do.call(rbind, lapply(names(dmats), function(tax_name){
  dmat <- dmats[[tax_name]]
  omat <- sapply(outliers[[tax_name]], function(i) {i | outliers[[tax_name]]})
  data.frame(distance=lowerTriangle(dmat), outlier=lowerTriangle(omat))
}))

dists$tax_name <- factor(rep(names(dmats), nseqs*(nseqs-1)/2))
with(dists, table(tax_name, outlier))


###################################################
### chunk number 13: 
###################################################
#line 195 "vignettes/clstutils/inst/doc/refSet.Rnw"
plot(bwplot(distance ~ tax_name, data=dists, ylim=c(0,0.15)))


###################################################
### chunk number 14: 
###################################################
#line 199 "vignettes/clstutils/inst/doc/refSet.Rnw"
plot(bwplot(distance ~ tax_name, data=subset(dists, !outlier), ylim=c(0,0.15)))


###################################################
### chunk number 15: 
###################################################
#line 209 "vignettes/clstutils/inst/doc/refSet.Rnw"
with(seqdat, {
  dmat <- ape::dist.dna(seqs, pairwise.deletion=TRUE, as.matrix=TRUE, model='raw')
  clstutils::prettyTree(nj(dmat), groups=tax_name,
                        ## type='unrooted',
                        X=outlier, fill=outlier)
})



###################################################
### chunk number 16: 
###################################################
#line 231 "vignettes/clstutils/inst/doc/refSet.Rnw"
with(seqdat[Efaecium,], {
  selected <- clstutils::maxDists(dmat, idx=which(isType),
                                  N=10, exclude=outlier, include.center=TRUE)
  prettyTree(nj(dmat), groups=ifelse(outlier,'outlier','non-outlier'),
             X=outlier,
             O=selected, fill=selected,
             labels=ifelse(isType,gettextf('type strain (%s)', accession),NA))
})