```##' neighbor-joining method
##'
##'
##' @title NJ
##' @param X distance matrix
##' @return phylo object
##' @author ygc
##' @examples
##' \dontrun{
##' X <- matrix(c(0,5,4,7,6,8,
##'		5,0,7,10,9,11,
##'		4,7,0,7,6,8,
##'		7,10,7,0,5,9,
##'		6,9,6,5,0,8,
##'		8,11,8,9,8,0), ncol=6)
##' rownames(X) <- colnames(X) <- LETTERS[1:6]
##' tree <- NJ(X)
##' print(tree)
##' }
NJ <- function(X) {
labels <- colnames(X)
N <- ncol(X)
otu_labs <- 1:N

dm <- as.matrix(X)
S <- colSums(dm)

## edge list of node 1 and node 2
edge1 <- edge2 <- numeric(2*N-3)
edge_length <- numeric(2*N-3)
k <- 1
cur_node <- 2*N-2
while (N > 3) {
ds <- 1e50
for (i in 1:(N-1)) {
for (j in (i+1):N) {
A <- N * dm[i,j] - S[i] - S[j]
if (A < ds) {
OUT1 <- i;
OUT2 <- j;
ds <- A
}
}
}
edge2[k] <- otu_labs[OUT1]
edge2[k+1] <- otu_labs[OUT2]
edge1[k] <- edge1[k+1] <- cur_node
dij <- dm[OUT1, OUT2]
B <- (S[OUT1]-S[OUT2]) / (N-2)
edge_length[k] <- (dij + B)/2
edge_length[k+1] <- (dij - B)/2

ij <- 1
new_dist <- numeric(N-2)
## d_kn <- 1/2 * (d_ik + d_jk - d_ij)
for (i in 1:N) {
if (i == OUT1 || i == OUT2) next
x <- dm[i, OUT1]
y <- dm[i, OUT2]
new_dist[ij] <- 1/2 * (x+y-dij)
ij <- ij + 1
}
## update data
dm <- dm[-c(OUT1, OUT2), -c(OUT1, OUT2)]
dm <- rbind(dm, new_dist)
dm <- cbind(dm, c(new_dist, 0))
otu_labs <- otu_labs[-c(OUT1, OUT2)]
otu_labs <- c(otu_labs, cur_node)
rownames(dm) <- otu_labs
colnames(dm) <- otu_labs
S <- colSums(dm)
cur_node <- cur_node-1
k <- k+2
N <- N - 1
}

n <- length(edge1)
edge1[(n-2):n] <- cur_node
edge2[(n-2):n] <- otu_labs
edge_length[n-2] <- (dm[1,2]+dm[1,3]-dm[2,3])/2
edge_length[n-1] <- (dm[2,1]+dm[2,3]-dm[1,3])/2
edge_length[n] <- (dm[3,1]+dm[3,2]-dm[1,2])/2
obj <- list(edge=cbind(as.numeric(edge1), as.numeric(edge2)),
edge.length=edge_length,
tip.label=labels, Nnode=length(labels)-2)
class(obj) <- "phylo"
return(obj)
}
```