.packageName <- "GeneNT"
BEST.kendall <- function(gene.name, Q, cormin)
{
    if (class(dat) != "matrix")
    dat <- as.matrix(dat)    
    qy <- which(names(dat[,0][1:nrow(dat)])==gene.name)
    p.list <- array(NA, nrow(dat))
    cor.list <- array(NA, nrow(dat))
    for (j in 1:nrow(dat))
    {
        g<- cor.test(dat[qy,], dat[j,], alternative = "two.sided", method = "kendall", exact = F)
        cor.list[j] <- cor(dat[qy,], dat[j,], method = "kendall")
        p.list[j] <- g$p.value
    }
    ###Calculate p-values for all pairs
    p.name1 <- rep(names(dat[,0][qy]),dim(dat)[1])
    p.name2 <- names(dat[,0][1:dim(dat)[1]])
    p.name <- cbind(p.name1, p.name2)
    colnames(p.name) <- c("gene1", "gene2")
    p.name <- data.frame(p.name)
    p.index1 <- rep(qy, dim(dat)[1])  
    p.index2 <- seq(1:dim(dat)[1])
    p.index <- cbind(p.index1, p.index2)
    colnames(p.index) <- c("index1", "index2")
    p.results <- cbind(p.index, p.name, cor.list, p.list)
    ###calculate q value
    fdr.out <- fdr.control(p.list, Q)
    q.list <- fdr.out$qvalues
    q.results <- cbind(p.results, q.list)
    ###use stepdown procedure to generate G1 set controlling FDR
    G1.results <- q.results[(q.results[,7] < Q),]
    sort.idx <- order(-abs(G1.results[, 5]))
    G1.results <- G1.results[sort.idx, ]
    G1 <- nrow(G1.results)
    pq <- G1/length(p.list)
    alpha = Q*pq
    ###calculate confidence intervals for G1
    lower <- array(NA, G1)
    higher <- array(NA, G1)
    
    if(G1 == 0)
    warning("Stage I screening returns no results. The Q value may be set to conservative (low)!")		
    for(i in 1:G1)
    {
        x <- G1.results[i,1] #get indices
        y <- G1.results[i,2]
        x <- dat[x,]
        y <- dat[y,]
        g <- kendall.confint(x,y,alpha)
        lower[i] <- g[1]
        if(g[2] >1)
            higher[i] <- 1
        if(g[2] <=1)
            higher[i] <- g[2]
    }
    lower <- as.matrix(lower)
    lower <- data.frame(lower)
    higher <- as.matrix(higher)
    higher <- data.frame(higher)    
    bkG1 <- cbind(G1.results, lower, higher) 
    indxg2 <- apply(bkG1,1,function(x) ifelse(((as.numeric(x[8]) > cormin) | (as.numeric(x[9]) < -cormin)), indxg2<- T, indxg2 <- F))
    bkG2 <- bkG1[indxg2,]
    write.table(bkG2, sep = "\t", file = "bkG2.txt")
    write.table(bkG1, sep = "\t", file = "bkG1.txt")
    cat("The screened pairs are now in your working directory. ")
    list(bkG1 = bkG1, bkG2 = bkG2)
}

BEST.pearson <- function(gene.name, Q, cormin, method = c("simple", "jackknife"))
{
    if (class(dat) != "matrix")
    dat <- as.matrix(dat)
    qy <- which(names(dat[,0][1:nrow(dat)])==gene.name)
    ##########
    method = match.arg(method)
    p.list <- array(NA, nrow(dat))
    cor.list <- array(NA, nrow(dat))
    if (method == "simple")
    {
        for (j in 1:nrow(dat))
        {
        g<- cor.test(dat[qy,], dat[j,], alternative = "two.sided", method = "pearson")
        cor.list[j] <- cor(dat[qy,], dat[j,])
        p.list[j] <- g$p.value
        }
    }
    if (method == "jackknife")
    {
        cor.temp <- array(NA, ncol(dat))
        for (j in 1:nrow(dat))
        {
            for (k in 1:ncol(dat))
            {
            if((sd(dat[qy,-k])*sd(dat[j,-k]))==0)
            cor.temp[k] <- 0
            if((sd(dat[qy,-k])*sd(dat[j,-k]))!=0)
            cor.temp[k] <- cor(dat[qy,-k], dat[j, -k], method = "pearson")
            }
        cor.list[j] <- median(cor.temp) 
        g<- cor.test(dat[qy,], dat[j,], alternative = "two.sided", method = "pearson")
        p.list[j] <- g$p.value
        }
    }
    ####calculate and rank p-values for all pairs
    p.name1 <- rep(names(dat[,0][qy]),dim(dat)[1])
    p.name2 <- names(dat[,0][1:dim(dat)[1]])
    p.name <- cbind(p.name1, p.name2)
    colnames(p.name) <- c("gene1", "gene2")
    p.name <- data.frame(p.name)
    p.index1 <- rep(qy, dim(dat)[1])  
    p.index2 <- seq(1:dim(dat)[1])
    p.index <- cbind(p.index1, p.index2)
    colnames(p.index) <- c("index1", "index2")
    p.results <- cbind(p.index, p.name, cor.list, p.list)

    ###calculate q values (FDR p-values)
    fdr.out <- fdr.control(p.list, Q)
    q.list <- fdr.out$qvalues
    q.results <- cbind(p.results, q.list)

    ###use stepdown procedure to generate G1 set controlling FDR
    G1.results <- q.results[(q.results[,7] < Q),]
    sort.idx <- order(-abs(G1.results[, 5]))
    G1.results <- G1.results[sort.idx, ]
    G1 <- nrow(G1.results)
    pq <- G1/length(p.list)

    ###construct confidence intervals for G1
    lower <- array(NA, G1)
    higher <- array(NA, G1)
    if(G1 == 0)
    warning("Stage I screening returns no results. The Q value may be set to conservative (low)!")		
    for(i in 1:G1)
    {
        x <- G1.results[i,1]
        y <- G1.results[i,2]
        x.row <- dat[x,]
        y.row <- dat[y,]
        alpha = Q*pq
        g <- cor.test(x.row, y.row, alternative = "two.sided", method = "pearson", conf.level = 1- alpha/2)
        lower[i] <- g$conf.int[1]
        higher[i] <- g$conf.int[2]
     }
    lower <- as.matrix(lower)
    lower <- data.frame(lower)
    higher <- as.matrix(higher)
    higher <- data.frame(higher)
    bpG1 <- cbind(G1.results, lower, higher) 
    indxg2 <- apply(bpG1,1,function(x) ifelse(((as.numeric(x[8]) > cormin) | (as.numeric(x[9]) < -cormin)), indxg2<- T, indxg2 <- F))
    bpG2 <- bpG1[indxg2,]
    write.table(bpG2, sep = "\t", file = "bpG2.txt")
    write.table(bpG1, sep = "\t", file = "bpG1.txt")
    cat("The screened pairs are now in your working directory. ")
    list(bpG1 = bpG1, bpG2 = bpG2)

}

cor.confint <- function (cor, N, alpha) 
{
        z <- atanh(cor)
        se <- 1/sqrt(N - 3)
        conf.int1 <- tanh(z-qnorm(1- alpha/2)*se)        
        conf.int2 <- tanh(z+qnorm(1- alpha/2)*se)        
	list(conf.int1 = conf.int1, conf.int2 = conf.int2)
}

corfdrci.inv <- function(cormin)
{
	#dat <- read.table(file.name, h = T, row.names = 1)
        if (class(dat) != "matrix")
    	dat <- as.matrix(dat)	
	gal.cor.m <- cor(t(dat))
	cov.p <- matrix(NA,nrow(dat), nrow(dat))

	for (i in 1:(nrow(dat)-1))
	{
		for (j in (i+1):nrow(dat))
		{
		g<- cor.test(dat[i,], dat[j,], alternative = "two.sided", method = "pearson")
		cov.p[i,j] <- g$p.value
		}
	}
	###Calculate p values for all pairs
	p.index <- sm.indexes(cov.p, diag = F)
	colnames(p.index) <- c("gene1", "gene2")
	cor.list <- sm2vec(t(gal.cor.m), diag = F) 
	p.list<- sm2vec(t(cov.p), diag = F)
	G.results <- cbind(p.index, p.list, cor.list)

	###step (2), sort the p values
	sort.index <- order(G.results[,3])
	G.results <- G.results[sort.index, ]
	###step (3),finding the min alpha  
	high <- array(NA, 99)
	low <- array(NA, 99)
	minalpha <- array(NA,length(cor.list))
	fdrp <- array(NA,length(cor.list))

	for (i in 1:length(cor.list))
	{
		x <- G.results[i,1]
		y <- G.results[i,2]

		#within cormin
		if(-cormin <= G.results[i,4] && G.results[i,4] <= cormin)
		minalpha[i] = 1
		
		#below neg cormin
		if(G.results[i,4] < -cormin)
		{
		for (al in 1:99)
		high[al] <- (cor.test(dat[x,], dat[y,], alternative = "two.sided", method = "pearson", conf.level = al/100))$conf.int[2] 
		
		if ((max(high) >= -cormin)&&(any(high < - cormin)))
		minalpha[i] <- min(which(high < -cormin))/100
		
		if(max(high) < -cormin)
		minalpha[i] <- 0
		
		if(!any(high < -cormin))
		minalpha[i] <- 1
		}

		#above pos cormin
		if(G.results[i,4] > cormin)
		{
		for (al in 1:99)
		low[al] <- (cor.test(dat[x,], dat[y,], alternative = "two.sided", method = "pearson", conf.level = al/100))$conf.int[1] 

		if((min(low) <= cormin)&&(any(low > cormin)))
		minalpha[i] <- (100- max(which(low > cormin)))/100
		
		if((min(low) > cormin))
		minalpha[i] <- 0

		if(!any(low > cormin))
		minalpha[i] <- 1
		}
	}  
	###step (4), compute the index
	boo <- array(NA, length(cor.list))
	fdrp <- array(NA, length(cor.list)) 
	for (j in 1:length(cor.list))
	{
		for (k in 1:length(cor.list))
		{
		boo[k] <- (G.results[k,3]*k/length(cor.list) <= minalpha[j])
		}
	fdrp[j] <- ifelse(sum(boo) == 0, G.results[1,3], G.results[sum(boo),3])
	}
	list(fdrp = fdrp)
}

corfdrci <- function(Q, cormin)
{
	if (class(dat) != "matrix")
    	dat <- as.matrix(dat)		
	gal.cor.m <- cor(t(dat))
	cov.p <- matrix(NA,nrow(dat), nrow(dat))
	for (i in 1:(nrow(dat)-1))
	{
		for (j in (i+1):nrow(dat))
		{
		g<- cor.test(dat[i,], dat[j,], alternative = "two.sided", method = "pearson")
		cov.p[i,j] <- g$p.value
		}
	}
	###Calculate p-values for all pairs
	p.name <- sm.name(gal.cor.m)
	colnames(p.name) <- c("gene1", "gene2")
	p.name <- data.frame(p.name)
	p.index <- sm.indexes(cov.p, diag = F)
	colnames(p.index) <- c("index1", "index2")
	p.list<- sm2vec(t(cov.p), diag = F)
	cor.list <- sm2vec(t(gal.cor.m), diag = F) 
	p.results <- cbind(p.index, p.name, cor.list, p.list)
	
	###calculate q values
	fdr.out <- fdr.control(p.list, Q)
	q.list <- fdr.out$qvalues
	q.results <- cbind(p.results, q.list)

	###use stepdown procedure to generate G1 set at certain FDR
	G1.results <- q.results[(q.results[,7] < Q),]
	sort.idx <- order(-abs(G1.results[, 5]))
	G1.results <- G1.results[sort.idx, ]
	G1 <- nrow(G1.results)
	pq <- G1/length(p.list)

	###calculate confidence intervals for G1
	lower <- array(NA, G1)
	higher <- array(NA, G1)
	alpha = Q*pq

	if(G1 == 0)
	warning("Stage I screening returns no results. The Q value may be set to conservative (low)!")		
	for(i in 1:G1)
	{
	x <- G1.results[i,1]
	y <- G1.results[i,2]
	r<- cor(dat[x,], dat[y,])
	g <- cor.confint(r, ncol(dat), alpha)
	lower[i] <- g$conf.int1
	higher[i] <- g$conf.int2
	}
	lower <- as.matrix(lower)
	lower <- data.frame(lower)
	higher <- as.matrix(higher)
	higher <- data.frame(higher)
	pG1 <- cbind(G1.results, lower, higher) 
	indxg2 <- apply(pG1,1,function(x) ifelse(((as.numeric(x[8]) > cormin) | (as.numeric(x[9]) < -cormin)), indxg2<- T, indxg2 <- F))
	pG2 <- pG1[indxg2,]
	write.table(pG2, sep = "\t", file = "pG2.txt")
	write.table(pG1, sep = "\t", file = "pG1.txt")
        cat("The screened pairs are now in your working directory. ")
	list(pG1=pG1, pG2=pG2)
}

getBM <- function(pG2, kG2)
{
   ##combine two screened pairs
   ppair <- pG2[,3:5]
   kpair <- kG2[,3:5]
   pair <- rbind(ppair, kpair)

   #find unsorted unique list of probsets
   index.us <- unique(rbind(as.matrix(pair[,1]), as.matrix(pair[,2])))
   #construct Boolean connectivity matrix to export to Pajek
   BM <- matrix(0, nrow(index.us), nrow(index.us)) 
   for(i in 1:nrow(pair))
   {
 	x <- which(index.us == as.character(pair[i,1])) 
 	y <- which(index.us == as.character(pair[i,2])) 
 	BM[x,y] <- 1
 	BM[y,x] <- 1
   }
   diag(BM) <- 1
   row.names(BM) <- index.us
   colnames(BM) <- index.us
   write.table(BM, sep = "\t", file = "BM.tsv")
   
   #the code belwo is obtained from http://vlado.fmf.uni-lj.si/pub/networks/pajek/howto/HowToR.htm
   #to save ordinary matrix in R to pajek compatible
   savematrix <- function(n,direct,twomode=1){
    if ((dim(n)[1] == dim(n)[2]) & (twomode!=2))
      { write(paste("*Vertices",dim(n)[1]), file = direct);
            write(paste(seq(1,length=dim(n)[1]),' "',rownames(n),
                  '"',sep=""), file = direct,append=TRUE);
            write("*Matrix", file = direct,append=TRUE);
            write(t(n),file = direct,ncolumns=dim(n)[1],
                  append=TRUE) }
    else
      { write(paste("*Vertices",sum(dim(n)),dim(n)[1]),
              file = direct);
            write(paste(1:dim(n)[1],' "',rownames(n),'"',sep=""),
                  file = direct,append=TRUE);
            write(paste(seq(dim(n)[1]+1,length=dim(n)[2]),' "',
                  colnames(n),'"',sep=""), file = direct,append=TRUE);
            write("*Matrix", file = direct, append=TRUE);
            write(t(n),file = direct, ncolumns=dim(n)[2],append=TRUE)}
      }   
   savematrix(BM,"BMPajek.mat")
}

kendall.confint <- function (x, y, alpha) 
{
	##the followed code is adopted in large part from http://www.stat.umn.edu/geyer/5601/examp/corr.html
        signs <- sign(outer(x, x, "-") * outer(y, y, "-"))
	tau <- mean(signs[lower.tri(signs)])
	#tau <- cor(x,y,method = "kendall")
	cvec <- apply(signs, 1, sum)
	n <- length(cvec)
	sigsq <- (2 / (n * (n - 1))) * (((2 * (n - 2)) / (n * (n - 1))) * var(cvec) + 1 - tau^2)
	zcrit <- qnorm(1 - alpha/2)
	conf.int <- tau + c(-1, 1) * zcrit * sqrt(sigsq)
	conf.int
}

kendallfdrci <- function(Q, cormin)
{
	if (class(dat) != "matrix")
    	dat <- as.matrix(dat)	
	gal.cor.m <- cor(t(dat), method = "kendall")
	cov.p <- matrix(NA,nrow(dat), nrow(dat))
	for (i in 1:(nrow(dat)-1))
	{
		for (j in (i+1):nrow(dat))
		{
		g <- cor.test(dat[i,], dat[j,], alternative = "two.sided", method = "kendall", exact = F)
		cov.p[i,j] <- g$p.value
		}
	}
	###list p-values for all pairs
	p.name <- sm.name(gal.cor.m)
	colnames(p.name) <- c("gene1", "gene2")
	p.name <- data.frame(p.name)
	p.index <- sm.indexes(cov.p, diag = F)
	colnames(p.index) <- c("index1", "index2")
	p.list<- sm2vec(t(cov.p), diag = F)
	cor.list <- sm2vec(t(gal.cor.m), diag = F) 
	p.results <- cbind(p.index, p.name, cor.list, p.list)

	###calculate q values 
	fdr.out <- fdr.control(p.list, Q)
	q.list <- fdr.out$qvalues
	q.results <- cbind(p.results, q.list)

	###use stepdown procedure to generate G1 set controlling FDR
	G1.results <- q.results[(q.results[,7] < Q),]
	sort.idx <- order(-abs(G1.results[, 5]))
	G1.results <- G1.results[sort.idx, ]
	G1 <- nrow(G1.results)
	pq <- G1/length(p.list)

	###calculate confidence intervals for G1
	lower <- array(NA, G1)
	higher <- array(NA, G1)
	alpha = Q*pq
	
	if(G1 == 0)
        warning("Stage I screening returns no results. The Q value may be set to conservative (low)!")		
	for(i in 1:G1)
	{
	x <- G1.results[i,1]
	y <- G1.results[i,2]
	x <- dat[x,]
	y <- dat[y,]
	g <- kendall.confint(x,y,alpha)
	lower[i] <- g[1]
	if(g[2] >1)
	higher[i] <- 1
	if(g[2] <=1)
	higher[i] <- g[2]
	}

	lower <- as.matrix(lower)
	lower <- data.frame(lower)
	higher <- as.matrix(higher)
	higher <- data.frame(higher)
	
	kG1 <- cbind(G1.results, lower, higher) 
	indxg2 <- apply(kG1,1,function(x) ifelse(((as.numeric(x[8]) > cormin) | (as.numeric(x[9]) < -cormin)), indxg2<- T, indxg2 <- F))
	kG2 <- kG1[indxg2,]
	write.table(kG2, sep = "\t", file = "kG2.txt")
	write.table(kG1, sep = "\t", file = "kG1.txt")
	cat("The screened pairs are now in your working directory. ")
	list(kG1=kG1, kG2=kG2) 
}

pcor.confint <- function (pcor, kappa, alpha) 
{
        z <- atanh(pcor)
        se <- 1/sqrt(kappa - 2)
        conf.int1 <- tanh(z-qnorm(1- alpha/2)*se)        
        conf.int2 <- tanh(z+qnorm(1- alpha/2)*se)        
	list(conf.int1 = conf.int1, conf.int2 = conf.int2)
}

pcorfdrci <- function(Q, pcormin)
{
	#dat <- read.table(file.name, h = T, row.names = 1)
        if (class(dat) != "matrix")
    	dat <- as.matrix(dat)	
	gal.cor.m <- cor(t(dat))
	b.pcor <- ggm.estimate.pcor(t(dat), method = "bagged.pcor")
	pcor.list <- sm2vec(b.pcor)
	
	pcov.p <- matrix(NA,nrow(dat), nrow(dat))
	kappa <- cor.fit.mixture(pcor.list)$kappa 
	for (i in 1:(nrow(dat)-1))
	{
		for (j in (i+1):nrow(dat))
		{
		r<- b.pcor[i,j]
		pval<- cor0.test(r, kappa, method = "ztransform")
		pcov.p[i,j] <- pval
		}
	}	
	###Calculate p values for all pairs
	p.name <- sm.name(gal.cor.m)
	colnames(p.name) <- c("gene1", "gene2")
	p.name <- data.frame(p.name)
	p.index <- sm.indexes(pcov.p, diag = F)
	colnames(p.index) <- c("index1", "index2")
	p.list<- sm2vec(t(pcov.p), diag = F)
	p.results <- cbind(p.index,p.name, pcor.list, p.list)

	###calculate q values
	fdr.out <- fdr.control(p.list, Q)
	q.list <- fdr.out$qvalues
	q.results <- cbind(p.results, q.list)

	###stepdown procedure to generate G1 subset controling FDR
	G1.results <- q.results[q.results[,7] < Q,]
	sort.idx <- order(-abs(G1.results[, 5]))
	G1.results <- G1.results[sort.idx, ]
	G1 <- nrow(G1.results)
	pq <- G1/length(p.list)

	###calculate confidence intervals for G1
	lower <- array(NA, G1)
	higher <- array(NA, G1)
	alpha = Q*pq

	for(i in 1:G1)
	{
		x <- G1.results[i,1]
		y <- G1.results[i,2]
		pcor <- b.pcor[x,y]
		g<- pcor.confint(pcor, kappa, alpha)

		lower[i] <- g$conf.int1
		higher[i] <- g$conf.int2
	}
	
	lower <- as.matrix(lower)
	lower <- data.frame(lower)
	higher <- as.matrix(higher)
	higher <- data.frame(higher)

	G1.all <- cbind(G1.results, lower, higher) 
	indxg2 <- apply(G1.all,1,function(x) ifelse(((as.numeric(x[8]) > pcormin) | (as.numeric(x[9]) < - pcormin)), indxg2<- T, indxg2 <- F))
	G2 <- G1.all[indxg2,]
	write.table(G2, sep = "\t", file = "G2.txt")
	write.table(G1.all, sep = "\t", file = "G1.txt")
	cat("The screened pairs are now in your working directory. ")
}

pcorfdrci.inv <- function(pcormin)
{
	#dat <- read.table(file.name, h = T, row.names = 1)
        if (class(dat) != "matrix")
    	dat <- as.matrix(dat)		
	###step (1), get the p values
	gal.cor.m <- cor(t(dat))
	b.cor <- bagged.cor(gal.cor.m)
	b.pcor <- cor2pcor(b.cor)
	pcor.list <- sm2vec(b.pcor)
	kappa <- cor0.estimate.kappa(pcor.list) 
	pcov.p <- matrix(NA,nrow(dat), nrow(dat))

	for (i in 1:(nrow(dat)-1))
	{
		for (j in (i+1):nrow(dat))
		{
		r<- b.pcor[i,j]
		pval<- cor0.test(r, kappa, method = "ztransform")
		pcov.p[i,j] <- pval
		}
	}
	###Calculate p values for all pairs
	p.index <- sm.indexes(pcov.p, diag = F)
	colnames(p.index) <- c("index1", "index2")
	p.list<- sm2vec(t(pcov.p), diag = F)
	G.results <- cbind(p.index, p.list, pcor.list)

	###step (2), sort the p values
	sort.index <- order(G.results[,3])
	G.results <- G.results[sort.index, ]
	
	###step (3),finding the min alpha  
	high <- array(NA, 99)
	low <- array(NA, 99)
	minalpha <- array(NA,length(pcor.list))
	fdrp <- array(NA,length(pcor.list))

	for (i in 1:length(pcor.list))
	{
	#within pcormins
	if(-pcormin <= G.results[i,4] && G.results[i,4] <= pcormin)
	minalpha[i] = 1

	#below neg pcormin
	if(G.results[i,4] < -pcormin)
	{
		for (al in 1:99)
		{
		g<- pcor.confint(G.results[i,4], kappa, al/100)
		high[al] <- g$conf.int2
		} 
		if ((max(high) >= -pcormin)&&(any(high < - pcormin)))
		minalpha[i] <- min(which(low > pcormin))/100

		if(max(high) < -pcormin)
		minalpha[i] <- 0

		if(!any(high < - pcormin))
		minalpha[i] <- 1
	}

	#above pos pcormin
	if(G.results[i,4] > pcormin)
	{
		for (al in 1:99)
		{
		g<- pcor.confint(G.results[i,4], kappa, al/100)
		low[al] <- g$conf.int1
		} 
		if((min(low) <= pcormin)&&(any(low > pcormin)))
		minalpha[i] <- min(which(low > pcormin))/100

		if((min(low) > pcormin))
		minalpha[i] <- 0

		if(!any(low > pcormin))
		minalpha[i] <- 1
	}
	}  
	###step (4), compute the index
	boo <- array(NA, length(pcor.list))
	fdrp <- array(NA, length(pcor.list)) 
	for (j in 1:length(pcor.list))
	{
		for (k in 1:length(pcor.list))
		{
		boo[k] <- (G.results[k,3]*k/length(pcor.list) <= minalpha[j])
		}
		fdrp[j] <- ifelse(sum(boo) == 0, G.results[1,3], G.results[sum(boo),3])
	}
        list(fdrp = fdrp)
}

priorclust <- function(p)
{
   #dat <- read.table(file.name, h = T,row.names = 1) 
   C <- cor(t(dat))
   diag(C) <- 0
   D <- (1-abs(C))^p
   diag(D) <- 0
   write.table(D, sep = "\t", file = "D.tsv")
   d <- as.dist(D)
   obj <- hclust(d)
   plot(obj, label = F, hang = 0, main = "Clustering with prior distance matrix", xlab = "" )
   y <- identify(obj, MAXCLUSTER = 50)
}

sm.name <- function (m) 
{
    l <- nrow(m)
    index1 <- rep(NA, l * (l - 1)/2)
    index2 <- rep(NA, l * (l - 1)/2)
    k <- 1
    for (i in 1:(l - 1)){ 
	for (j in (i + 1):l) {
        index1[k] <- names(m[i])
        index2[k] <- names(m[j])
        k <- k + 1
    }
    }
    return(cbind(index1, index2))
}

spclust <- function(p, pG2, kG2)
{
   ##combine two screened pairs
   ppair <- pG2[,3:5]
   kpair <- kG2[,3:5]
   pair <- rbind(ppair, kpair)

   #find unsorted unique list of probsets
   index.us <- unique(rbind(as.matrix(pair[,1]), as.matrix(pair[,2])))

   #calculate the connectivity list
   mdegree <- matrix(0, nrow(index.us), nrow(index.us)) 
   for(i in 1:nrow(pair))
   {
 	x <- which(index.us == as.character(pair[i,1])) 
 	y <- which(index.us == as.character(pair[i,2])) 
 	mdegree[x,y] <- 1
 	mdegree[y,x] <- 1
   }

   #order probsets according to connecitvity
   ll <- apply(mdegree, 1, sum)
   ll <- as.matrix(ll)
   degreedata <- cbind(index.us, ll)
   idx <- order(as.numeric(degreedata[,2]), decreasing = T)
   index <- degreedata[idx,1]
   index <- data.frame(index)

   M <- matrix(NA, nrow(index), nrow(index)) 
   for(i in 1:nrow(pair))
   {
 	x <- which(index == as.character(pair[i,1])) 
 	y <- which(index == as.character(pair[i,2])) 
        if(is.na(M[x,y])) 
  	{
   	   M[x,y] <- (1- abs(pair[i,3]))^p
           M[y,x] <- (1- abs(pair[i,3]))^p
 	} 
 	if(!is.na(M[x,y])) #choose the largest correlation, shortest path
 	{ 
   	   M[x,y] <- min((1- abs(pair[i,3]))^p, M[x,y])   
   	   M[y,x] <- M[x,y]
 	}     
   }

   rcname <- as.matrix(index)
   colnames(M) <- rcname
   row.names(M) <- rcname
   diag(M) <- 0

   z1 <- allShortestPaths(M)
   #mdegree is the distance matrix filled with SPs (SP distance matrix)
   Mdist <- z1$length
   colnames(Mdist) <- rcname
   row.names(Mdist) <- rcname

   #find the GSC.
   if(!any(is.na(Mdist)))
   { gsc <- Mdist }
   if(any(is.na(Mdist)))
   { 
      i = 1
      while(1)
      {
 	 if(any(is.na(Mdist[1:i, 1:i])))
 	 break
 	 i <- i + 1 
      }
      gsc <- Mdist[1:(i-1), 1:(i-1)]
   }
   
   if(any(is.na(gsc))) #check whether NA is mising  
   cat("GCC search error!")
   
   write.table(gsc, sep = "\t", file = "gsc.tsv")
   d <- as.dist(gsc)
   obj <- hclust(d, method = "single")
   write.table(obj$label[obj$order], sep = "\t", file = "label.tsv")   
   cat("The Shortest-Path distance matrix and labels have been written to the default directory!")   
   plot(obj, label = F, hang = 0, main = "Clustering with posterior (shortest-path) distance matrix", sub = "", xlab = "" )
}

