---
output:
  html_document: default
  word_document: default
  pdf_document: default
---

```{r message=FALSE, warning=FALSE}
# library(openxlsx)
library(magrittr)
```

# KC Correlations Calculation 
```{r}
#define gauss.k kernel
gaussk <- function(u,v, sig){ exp(-sig*sum( (u-v)^2) ) }


# kc with centralizing the kernel matrix
corr.gene.gaussk <- function(timeT, matG, sig) {
  ### --------------------- ###
  # matG: 
  #   nrow: the number of features(genes); row.names of the matrix should be the gene name list
  #   ncol: the number of time points
  #
  # sig: 
  #   The parameter for RBF kernel function
  #
  ### --------------------- ###
  
  tlength <- length(timeT)
  
  # calculate kernel matrix K
  matK <- matrix(0, nrow = tlength, ncol = tlength)
  for (i in 1:tlength) {
    for (j in 1:tlength)
      matK[i, j] <- gaussk(matG[, i], matG[, j], sig)
  }
  
  # center kernel matrix
  IminusL <-
    diag(1, nrow = tlength, ncol = tlength) - (1 / tlength) * matrix(1, nrow =
                                                                       tlength, ncol = tlength)
  matK.center <- IminusL %*% matK %*% IminusL
  
  #calculate the dot product between matG and matK.center
  matGK <- matG %*%  matK.center
  
  #correlation
  return(cor(matGK[1, ], matGK[2, ]))
  
}

# kc calculation with centralizing the kernel matrix -- select the gene pair by users
kc_genes = function(timeT, mat, sig, gene_names = NULL) {
  ### --------------------- ###
  # matG:
  #   nrow: the number of features(genes); row.names of the matrix should be the gene name list
  #   ncol: the number of time points
  #
  # sig:
  #   The parameter for RBF kernel function
  #
  # gene_names:
  #   Return the kc of the specific gene pair(If `NULL`, return all results)
  ### --------------------- ###
  p = nrow(mat)
  
  if (is.null(gene_names)) {
    # gene_name_list = row.names(mat)
    # mat_kcs = matrix(nrow = p, ncol = p)
    ### --------------------- ###
    # matG:
    #   nrow: the number of features(genes); row.names of the matrix should be the gene name list
    #   ncol: the number of time points
    #
    # sig:
    #   The parameter for RBF kernel function
    #
    ### --------------------- ###
    matG = mat
    tlength <- length(timeT)
    plength = nrow(matG)
    mat_kc = matrix(nrow=plength, ncol=plength)
    
    for(u in 1:plength) {
      for (v in 1:plength) {
        
        # calculate kernel matrix K
        matK <- matrix(0, nrow = tlength, ncol = tlength)
        for (i in 1:tlength) {
          for (j in 1:tlength)
            matK[i, j] <- gaussk(matG[u, i], matG[v, j], sig)
        }
        
        # center kernel matrix
        IminusL <-
          diag(1, nrow = tlength, ncol = tlength) - (1 / tlength) * matrix(1, nrow = tlength, ncol = tlength)
        matK.center <- IminusL %*% matK %*% IminusL
        
        # Calculate the dot product between matG and matK.center
        matGK <- matG %*% matK.center
        # Correlations
        mat_kc[u, v] = (cor(matGK[1,], matGK[2,])) 
      }
    }
    
    row.names(mat_kc) = colnames(mat_kc) = row.names(matG)
    print("--- kc correlations calculation ---")
    return(mat_kc)
  }
  else {
    
    gene_name_list = row.names(mat)
    ids = which(gene_name_list %in% gene_names)
    genes = gene_name_list[ids]
    # print(paste0("The gene pair: ", genes))
    
    matG = mat[ids, ]
    tlength <- length(timeT)
    
    # calculate kernel matrix K
    matK <- matrix(0, nrow = tlength, ncol = tlength)
    for (i in 1:tlength) {
      for (j in 1:tlength)
        matK[i, j] <- gaussk(matG[, i], matG[, j], sig)
    }
    
    # center kernel matrix
    IminusL <-
      diag(1, nrow = tlength, ncol = tlength) - (1 / tlength) * matrix(1, nrow = tlength, ncol = tlength)
    matK.center <- IminusL %*% matK %*% IminusL
    
    # Calculate the dot product between matG & matK.center
    matGK <- matG %*% matK.center
    
    #correlation
    cor_kc = cor(matGK[1, ], matGK[2, ])
    # print(matGK)
    print(paste0("The kc correlation for the gene pair(", gene_names[1], ", ", gene_names[2], "): ", cor_kc))
    return(cor_kc)
    
  }
  
}

```

## Example (application 1, replication 1).
```{r}
# file_names = dir("Data/readcount/")[grep(".csv", dir("Data/readcount/"))]
# expr_data.DEseq = lapply(1:length(file_names), function(i) {
#     data = strsplit(readLines(paste0("Data/readcount/", file_names[i])), ",")
#     data = data[33:40]
#     return(t(matrix(unlist(data),  nrow = 9)))
# })
```

```{r}
File_name = "Data/RNA-seq__8_genes_data_ DESeq.xlsx"
Sheet_names = openxlsx::getSheetNames("Data/RNA-seq__8_genes_data_ DESeq.xlsx")[-1]
# Sheet names: "Th17_R1", "Th17_R2", "Th17_R3", "Th0_R1", "Th0_R2", "Th0_R3" 
expr_data.DEseq_ = lapply(1:length(Sheet_names), function(i) {
    data = openxlsx::read.xlsx(File_name, sheet = Sheet_names[i])
    data = data[16:23, 5:9]
    data = data %>% dplyr::mutate_all(function(x) as.numeric(x))
    return(data)
})
```

```{r}
DESeq_th0_R1 = expr_data.DEseq_[[4]]
# DESeq_th0_R1 = DESeq_th0_R1[, 5:9]
# DESeq_th0_R1 = DESeq_th0_R1 %>% as.numeric() %>%  matrix(nrow = 8, ncol = 5, byrow = F)
```

* Standardization
  + row = 8(genes), col = 5(timeT)
```{r}
gene_names = c("TIAM1", "RORC", "KIF11", "MAP1B", "ISG20", "IL17A", "IL17F", "RAB13")
DESeq_th0_R1.std = DESeq_th0_R1 %>% apply(MARGIN = 1, scale) %>% t() %>% set_rownames(gene_names)

print(kc_genes(1:5, DESeq_th0_R1.std, sig = 0.5, gene_names = c("TIAM1", "IL17A")))
```

