h2zconfirmation.Rmd

---
title: "h2zconfirmation"
output: html_document
date: "2024-05-19"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## libraries
```{r}
library(dendextend)
library(tidyverse)
library(dplyr)
```

# load and setup
```{r, results='hold', warning=FALSE, message=FALSE}
#zz sem lines
zzsl<- read_tsv("z2z-sem-lines.mat")
zzsl$...1 = str_replace_all(zzsl$...1, "_",":")
colnames(zzsl) = str_replace_all(colnames(zzsl), "_",":")
zgenes <- as.matrix(zzsl[,1])
zgenes <- str_replace_all(zgenes, "_",":")
data<-as.dist(zzsl[,-1])

#orthologs
orths = read.table("orthologs.txt", col.names = c("human", "z1","z2","z3","z4"), header = FALSE, sep = "\t", fill = TRUE)

#zz txt lines
zztl<- read_tsv("z2z-txt-lines.mat")
zztl$...1 = str_replace_all(zztl$...1, "_",":")
colnames(zztl) = str_replace_all(colnames(zztl), "_",":")
zgenestxt <- as.matrix(zzsl[,1])
zgenestxt <- str_replace_all(zgenestxt, "_",":")
datatxt<-as.dist(zztl[,-1])

#h2z sem lines
hzsl<- read_tsv("bigtable.sem-lines.tsv")
hzsl$zfgene = str_replace_all(hzsl$zfgene, "_",":")
colnames(hzsl) = str_replace_all(colnames(hzsl), "_",":")
hgenes <- colnames(hzsl)
zgeneshz <- as.matrix(hzsl[,1])
zgeneshz <- str_replace_all(zgeneshz, "_",":")

# h2z txtlines
hztl<- read_tsv("bigtable.txt-lines.tsv")
hztl$zfgene = str_replace_all(hztl$zfgene, "_",":")
colnames(hztl) = str_replace_all(colnames(hztl), "_",":")
zgeneshz <- as.matrix(hzsl[,1])
zgeneshz <- str_replace_all(zgeneshz, "_",":")



```

# cluster analysis z2z
```{r}
library(cluster)
set.seed(125)

myclustertxt <- function(x, k) list(cluster=cutree(hclust(datatxt, method = "average"),k=k))

myclusGaptxt <- clusGap(as.matrix(zztl[,-1]),
                     FUN = myclustertxt,
                     K.max = 350,
                     B = 2)

plot(myclusGaptxt, main = "Gap Statistic Text Comparison")+
abline(v=314, col = "red")

txtclustnum = with(myclusGaptxt, maxSE(Tab[,"gap"], Tab[,"SE.sim"], method="globalSEmax")) # 314



myclustersem <- function(x, k) list(cluster=cutree(hclust(data, method = "average"),k=k))

myclusGapsem <- clusGap(as.matrix(zzsl[,-1]),
                     FUN = myclustersem,
                     K.max = 350,
                     B = 2)

plot(myclusGapsem, main = "Gap Statistic Semantic Comparison")+
abline(v=348, col = "red")

semclustnum = with(myclusGapsem, maxSE(Tab[,"gap"], Tab[,"SE.sim"], method="globalSEmax")) # 348

```




#zztl
```{r}
hclusttxt <- hclust(datatxt, method = 'average')

cut1txt <- cutree(hclusttxt, h=.1)

cluster1txt <- zztl[,-1]
cluster1txt <- mutate(cluster1txt, cluster = cut1txt)
#head(cluster1)
counttxt <- count(cluster1txt, cluster)
summary(counttxt[2])
head(counttxt)

clusterstxt = as.data.frame(as.matrix(cut1txt))
clusterstxt$gene = zgenes
colnames(clusterstxt)[1] <- "cluster"
colnames(clusterstxt)[2] <- "gene"

head(clusterstxt)

dendtxt <- as.dendrogram(hclusttxt)
# k= no of clusters, h= height of cut
dendtxt <- color_branches(dendtxt, h=0.1)
plot(dendtxt)
plot(dendtxt, ylab = "Semantic similarity score", main = "Dendogram of Zebrafish to Zebrafish Text Comparison Scores")




```



# plot 0.1
```{r}
hclust <- hclust(data, method = 'average')
dend <- as.dendrogram(hclust)
# k= no of clusters, h= height of cut
dend <- color_branches(dend, h=0.1)
plot(dend, ylab = "Semantic similarity score", main = "Dendogram of Zebrafish to Zebrafish Semantic Scores")

```
# cuts 0.1
```{r}
cut1 <- cutree(hclust, h=.1)

cluster1 <- zzsl[,-1]
cluster1 <- mutate(cluster1, cluster = cut1)
#head(cluster1)
count <- count(cluster1, cluster)
summary(count[2])
head(count)

```
# add gene names to cluster matrix -- clusters
```{r}
clusters = as.data.frame(as.matrix(cut1))
clusters$gene = zgenes
colnames(clusters)[1] <- "cluster"
colnames(clusters)[2] <- "gene"

head(clusters)

```


# pull out
```{r}
# target = "SLX4"
# score = matrix()
# 
# 
# b = hzsl %>%
#   select(T = target,  zfgene)
# 
# tops = arrange(b, T)[1:5,]
# tops
# 
# #match=orths[which(orths == target, arr.ind = TRUE)[1],1]
# 
# 
# 
# for(i in 1:5) {
#   i = tops[[i,2]]
#   if (orths[which(orths == i, arr.ind = TRUE)[1],1] == target)
#     score[target] = 1;
# }
# 
# clustnum = clusters[which(clusters == i, arr.ind = TRUE)[1],1]
# clustnum
# subset = filter(clusters, cluster == clustnum)
# for(poss in subset[,2]){
#   if (orths[which(orths == poss, arr.ind = TRUE)[1],1] == target)
#     score[2]=1
# }
# 
# score
```


```{r}


```

# auto
```{r}
# score matriz - for each h gene, ranking of ortholog and whether ortholog in cluster of highest rank
score = matrix(0, nrow = 647, ncol=7, dimnames = list(hgenes[-1],c("gene", "rankingSem","captureSem.1","captureSemOpt", "rankingTxt","captureTxt.1", "captureTxtOpt")))
score[,1]= hgenes[-1]

semorthscores = vector()


# for each zgene sem
for (htarget in hgenes[-1]) {
  ###### change to h2z #############
  
  # seperate out 8 smallest scores and name of associated gene
  b = hzsl %>%
  select(t=htarget, zfgene)
  tops = arrange(b, t)
  zmatches = orths[(which(orths == htarget)[1]),-1]
  #target = str_replace(target, "_",":")
  
  
  
  # for the top 8
  for(znum in 1:731) {
    
  # j is the gene in top 8
  zname = tops[[znum,2]]
  zname = str_replace(zname, "_",":")
  #print(j)
  
  # find ortholog of top scoring gene
  currscore = orths[which(orths == zname, arr.ind = TRUE)[1],1]
  #print(currscore)

  # if ortholog of top score is the same as ortholog of source gene, correctly identified match
  if (currscore == htarget){
    #print(target)
    score[htarget,2] = znum
    semorthscores = append(semorthscores,tops[znum,1])
  }
  

  if (score[htarget,3] == 0){
    #print(znum)
    # identify cluster of top ranking gene
    clustnum = clusters[which(clusters == zname, arr.ind = TRUE)[1],1]
    # list all genes in that cluster
    subset = filter(clusters, cluster == clustnum)
    # if cluster of top ranking gene includes the source gene, correctly captured
    for(poss in subset[,2]){
      if(poss %in% zmatches)
        score[htarget,3] = znum

    }
  }
  
  
}
}
```


```{r}

txtorthscores = vector()

for (htarget in hgenes[-1]) {
  ###### change to h2z #############
  
  # seperate out 8 smallest scores and name of associated gene
  b = hztl %>%
  select(T = htarget,  zfgene)
  tops = arrange(b, T)
  zmatches = orths[(which(orths == htarget)[1]),-1]
  #target = str_replace(target, "_",":")
  
  
  # for the top 8
  for(znum in 1:731) {
    
  # j is the gene in top 8
  zname = tops[[znum,2]]
  zname = str_replace(zname, "_",":")
  #print(j)
  
  # find ortholog of top scoring gene
  currscore = orths[which(orths == zname, arr.ind = TRUE)[1],1]
  #print(currscore)
  
  txtorthscores = append(semorthscores,tops[1,1])

  # if ortholog of top score is the same as ortholog of source gene, correctly identified match
  if (currscore == htarget){
    #print(target)
    score[htarget,5] = znum
  }
  

  if (score[htarget,6] == 0){
    #print(znum)
    # identify cluster of top ranking gene
    clustnum = clusterstxt[which(clusterstxt == zname, arr.ind = TRUE)[1],1]
    # list all genes in that cluster
    subset = filter(clusterstxt, cluster == clustnum)
    # if cluster of top ranking gene includes the source gene, correctly captured
    for(poss in subset[,2]){
      if(poss %in% zmatches)
        score[htarget,6] = znum

    }
  }
  
  
}
}

S = data.frame(score)

```

# opt cluster
```{r}
cuttxtopt <- cutree(hclusttxt, k= 314)
clustertxtopt = as.data.frame(as.matrix(cuttxtopt))
clustertxtopt$gene = zgenes
colnames(clustertxtopt)[1] <- "cluster"
colnames(clustertxtopt)[2] <- "gene"

cutsemopt <- cutree(hclust, k= 348)
clustersemopt = as.data.frame(as.matrix(cutsemopt))
clustersemopt$gene = zgenes
colnames(clustersemopt)[1] <- "cluster"
colnames(clustersemopt)[2] <- "gene"


```

```{r}
for (htarget in hgenes[-1]) {
  ###### change to h2z #############
  
  # seperate out 8 smallest scores and name of associated gene
  b = hztl %>%
  select(T = htarget,  zfgene)
  tops = arrange(b, T)
  zmatches = orths[(which(orths == htarget)[1]),-1]
  #target = str_replace(target, "_",":")
  
  
  
  # for the top 8
  for(znum in 1:731) {
    
  # j is the gene in top 8
  zname = tops[[znum,2]]
  zname = str_replace(zname, "_",":")
  #print(j)

  if (score[htarget,7] == 0){
    #print(znum)
    # identify cluster of top ranking gene
    clustnum = clustertxtopt[which(clustertxtopt == zname, arr.ind = TRUE)[1],1]
    # list all genes in that cluster
    subset = filter(clustertxtopt, cluster == clustnum)
    # if cluster of top ranking gene includes the source gene, correctly captured
    for(poss in subset[,2]){
      if(poss %in% zmatches)
        score[htarget,7] = znum

    }
  }
  
  
  }
  
  # seperate out 8 smallest scores and name of associated gene
  b = hzsl %>%
  select(T = htarget,  zfgene)
  tops = arrange(b, T)
  zmatches = orths[(which(orths == htarget)[1]),-1]
  #target = str_replace(target, "_",":")
  
  
  # for the top 8
  for(znum in 1:731) {
    
  # j is the gene in top 8
  zname = tops[[znum,2]]
  zname = str_replace(zname, "_",":")
  #print(j)

  if (score[htarget,4] == 0){
    #print(znum)
    # identify cluster of top ranking gene
    clustnum = clustersemopt[which(clustersemopt == zname, arr.ind = TRUE)[1],1]
    # list all genes in that cluster
    subset = filter(clustersemopt, cluster == clustnum)
    # if cluster of top ranking gene includes the source gene, correctly captured
    for(poss in subset[,2]){
      if(poss %in% zmatches)
        score[htarget,4] = znum

    }
  }
  
  
  }
  
  
  
}
```


#summarize

```{r}
S %>% group_by(rankingSem) %>% summarise(ct = n())

ggplot(S, aes(x=rankingSem))+
  geom_bar(stat="count")
  #scale_x_discrete(limits=c("0","50","100","150","200","250","300","350","400","450","500","550","600"))

S %>% group_by(captureSem) %>% summarise(ct = n())

ggplot(S, aes(x=captureSem))+
  geom_histogram(stat = "count", binwidth = 1)

S %>% group_by(rankingTxt) %>% summarise(ct = n())
S %>% group_by(captureTxt) %>% summarise(ct = n())

#ggplot(S, aes(x=cluster))+
 # geom_histogram(stat = "count")
```

```{r}
S = data.frame(score)

S$rankingSem = as.numeric(S$rankingSem)
S$captureSem = as.numeric(S$captureSem)
S$rankingTxt = as.numeric(S$rankingTxt)
S$captureTxt = as.numeric(S$captureTxt)

S$captureSemOpt = as.numeric(S$captureSemOpt)
S$captureTxtOpt = as.numeric(S$captureTxtOpt)

hist(S$rankingSem)
hist(S$captureSem)
hist(S$rankingTxt)
hist(S$captureTxt)

hist(S$captureSemOpt)
hist(S$captureTxtOpt)
```



# orth aves
```{r}
mean(as.data.frame(semorthscores)[,2])
mean(as.data.frame(txtorthscores)[,2])

```