library("Biostrings")
library(ape)
library(K2)

K2Dist <- function (data, q=2)
{
  K2s=matrix(0,ncol=data$total,nrow=data$total)
  
  for (i in 1:(data$total))
  {
    for (j in i:data$total)
    {
      K2s[i,j]=K2similarity(data$seq[i], data$seq[j], q)
      K2s[j,i]=K2s[i,j]
    }
  }
  
  colnames(K2s)=data$names
  rownames(K2s)=data$names
  K2s
}


K2sDist <- function (data, SigmaSize=4)
{
  K2s=matrix(0,ncol=data$total,nrow=data$total)
  
  for (i in 1:(data$total-1))
  {
    for (j in i:data$total)
    {
      K2s[i,j]=K2star(data$seq[i], data$seq[j], SigmaSize)
      K2s[j,i]=K2s[i,j]
    }
  }
  
  colnames(K2s)=data$names
  rownames(K2s)=data$names
  K2s
}



CRMexperiment <- function (read.seq)
{
  index=1
  TOP=c(300, 136, 253, 36, 136, 36, 300)
  for (f in read.seq$V1)
  {
    file=paste("./Data/",f, sep="")
    fasta=readFasta(file)
    
    for (q in 2:6)
    {
      K2=rep(0,fasta$total^2)	#Kendall Similarity
      sign=as.vector(c(rep(1,fasta$total/2),rep(-1,fasta$total/2)) %*% t(c(rep(1,fasta$total/2),rep(-1,fasta$total/2))))
      self=seq(1,fasta$total+1,fasta$total^2)
      
      count=1
      for (i in 1:fasta$total)
      {
        for (j in 1:fasta$total)
        {
          K2[count]=K2similarity(fasta$seq[i], fasta$seq[j], q)	
          count=count+1
        }
      }
      K2=K2[-self]
      sign=sign[-self]
      out=data.frame(rank(-K2,  ties.method = "first"),sign)
      names(out)=c("rank","sign")
      
      top=TOP[index]
      cat("K2","(q=",q,")",f,sum(out[out$rank<=top,2]==1)/length(out[out$rank<=top,2]),"\n")
    }
    
    K2=rep(0,fasta$total^2)	#Kendall Similarity
    sign=as.vector(c(rep(1,fasta$total/2),rep(-1,fasta$total/2)) %*% t(c(rep(1,fasta$total/2),rep(-1,fasta$total/2))))
    self=seq(1,fasta$total+1,fasta$total^2)
    
    count=1
    for (i in 1:fasta$total)
    {
      for (j in 1:fasta$total)
      {
        K2[count]=K2star(fasta$seq[i], fasta$seq[j], 4)	
        count=count+1
      }
    }
    K2=K2[-self]
    sign=sign[-self]
    out=data.frame(rank(-K2,  ties.method = "first"),sign)
    names(out)=c("rank","sign")
    
    top=TOP[index]
    cat("K2*",f,sum(out[out$rank<=top,2]==1)/length(out[out$rank<=top,2]),"\n")
    
    
    index=index+1
  }
}

readFas<-function (file) 
{
  lines = read.table(file, stringsAsFactors = F,sep="\t")
  len = nrow(lines)
  ind = rep(0, len)
  name = rep("", len)
  count = 1
  for (i in 1:len) {
    if (substr(lines[i, 1], 1, 1) == ">") {
      name[count] = substring(lines[i, 1], 2)
      ind[count] = i
      count = count + 1
    }
  }
  ind[count] = len + 1
  seq = rep("", count - 1)
  for (i in 2:count) {
    for (j in (ind[i - 1] + 1):(ind[i] - 1)) {
      seq[i - 1] = paste(seq[i - 1], lines[j, 1], sep = "")
    }
  }
  list(seq = seq, total = count - 1, names = name[1:(count - 1)])
}

DNAexperiment <- function(DNA.data, pdffile="DNA.pdf")
{
  cao = read.tree(text = "((platypus,(opossum,wallaroo)),((mouse,rat),((((fin_whale,bluewhale),cow ),((whiterhinoceros,horse),(cat,(grayseal,harborseal)))),(gibbon,((orangutan,orangutan_abelii),(gorilla,((baboon,common_chimpanzee),human)))))));")
  

  pdf(pdffile)
  for (qq in 2:9)
  {
    K2mat = K2Dist(DNA.data, q = qq)
    plot(
      hclust(as.dist(-K2mat)),
      main = paste("Method=", "K2", ",Q=", qq),
      sub = "",
      yaxt = "n",
      xlab = "",
      ylab = ""
    )
    row.names(K2mat) = DNAlist
    dist = as.dist(1 - K2mat)
    phy1 = as.phylo(hclust(dist))
    cat("K2 (q=", qq, ") RF distance=", dist.topo(cao, phy1), "\n")
  }
  
  K2mat = K2sDist(DNA.data, 4)
  plot(hclust(as.dist(-K2mat)),main=paste("Method=","K2*"),sub="",yaxt="n",xlab = "",  ylab = "")
  row.names(K2mat) = DNAlist
  dist = as.dist(1 - K2mat)
  phy1 = as.phylo(hclust(dist))
  cat("K2* RF distance=", dist.topo(cao, phy1), "\n")
  dev.off()
}


FishExperiment <- function(dat, pdffile="Fish.pdf")
{
  fishtree = read.tree(text = "((Alepocephalus_agassizii,Bajacalifornia_megalops),(((Parajulis_poecilepterus,Halichoeres_melanurus),(Pteragogus_flagellifer,(Pseudolabrus_eoethinus,Pseudolabrus_sieboldi))),(((Ditrema_temminckii,Cymatogaster_aggregata),(Abudefduf_vaigiensis,Amphiprion_ocellaris)),((Etroplus_maculatus,Paretroplus_maculatus),((Paratilapia_polleni,Ptychochromoides_katria),((Hypselecara_temporalis,Astronotus_ocellatus),(Tylochromis_polylepis,((Oreochromis_sp._KM-2006,(Oreochromis_niloticus,Oreochromis_aureus)),(Neolamprologus_brichardi,Tropheus_duboisi)))))))));")
  
  FishList = read.csv(
    "./Data/Fish.list",
    stringsAsFactors = F,
    header = F,
    sep = ","
  )[, 2]
  
  pdf(pdffile)
  for (qq in 2:9)
  {
    K2mat = K2Dist(dat, q = qq)
    row.names(K2mat) = FishList
    plot(
      hclust(as.dist(-K2mat)),
      main = paste("Method=", "K2", ",Q=", qq),
      sub = "",
      yaxt = "n",
      xlab = "",
      ylab = ""
    )
    k2s = as.phylo(hclust(as.dist(1 - K2mat)))
    
    cat("K2 (q=", qq, ") RF distance=", dist.topo(fishtree, k2s), "\n")
  }
  
  
  K2mat = K2sDist(dat, 4)
  row.names(K2mat) = FishList
  plot(hclust(as.dist(-K2mat)),main=paste("Method=","K2*"),sub="",yaxt="n",xlab = "",  ylab = "")
  k2s = as.phylo(hclust(as.dist(1 - K2mat)))
  cat("K2* RF distance=", dist.topo(fishtree, k2s), "\n")
  
  dev.off()
}

# DNA.data= readFasta("./Data/CaoDNA.txt")
# DNAexperiment(DNA.data)
# 
# 
# dat=readFas("./Data/FishData.fas")
# FishExperiment(dat)
# 
# seq=read.table("./Data/CRMs.lst",stringsAsFactors = F)
# CRMexperiment(seq)

