geneCoexAnalysisE17.5Yuzwa.qmd

---
title: "Gene Correlation Analysis E17.5"
---

## Prologue

```{r}
library(COTAN)
library(ComplexHeatmap)
library(circlize)
library(dplyr)
library(Hmisc)
library(Seurat)
library(patchwork)
library(Rfast)
library(parallel)
library(doParallel)
library(HiClimR)
library(stringr)
library(fst)

options(parallelly.fork.enable = TRUE)

dataSetFile <- ("Data/Yuzwa_MouseCortex/CorticalCells_GSM2861514_E175.cotan.RDS")

name <- str_split(dataSetFile,pattern = "/",simplify = T)[3]
name <- str_remove(name,pattern = ".RDS")

project = "E17.5"


setLoggingLevel(1)
outDir <- "CoexData/"
setLoggingFile(paste0(outDir, "Logs/",name,".log"))


obj <- readRDS(dataSetFile)
file_code = getMetadataElement(obj, datasetTags()[["cond"]])

```

```{r}
source("src/Functions.R")
```

To compare the ability of COTAN to asses the real correlation between genes we define some pools of genes:

1.  Constitutive genes
2.  Neural progenitor genes
3.  Pan neuronal genes
4.  Some layer marker genes

```{r}
genesList <- list(
  "NPGs"= 
    c("Nes", "Vim", "Sox2", "Sox1", "Notch1", "Hes1", "Hes5", "Pax6"),
  "PNGs"= 
    c("Map2", "Tubb3", "Neurod1", "Nefm", "Nefl", "Dcx", "Tbr1"),
  "hk"= 
    c("Calm1", "Cox6b1", "Ppia", "Rpl18", "Cox7c", "Erh", "H3f3a",
      "Taf1", "Taf2", "Gapdh", "Actb", "Golph3", "Zfr", "Sub1",
      "Tars", "Amacr"),
  "layers" = 
    c("Reln","Lhx5","Cux1","Satb2","Tle1","Mef2c","Rorb","Sox5","Bcl11b","Fezf2","Foxp2","Ntf3","Rasgrf2","Pvrl3", "Cux2","Slc17a6", "Sema3c","Thsd7a", "Sulf2", "Kcnk2","Grik3", "Etv1", "Tle4", "Tmem200a", "Glra2", "Etv1","Htr1f", "Sulf1","Rxfp1", "Syt6") 
  # From https://www.science.org/doi/10.1126/science.aam8999
)
```

## COTAN

```{r}
int.genes <-getGenes(obj)
```

```{r}
coexMat.big <- getGenesCoex(obj)#[int.genes,int.genes]

coexMat <- getGenesCoex(obj)[c(genesList$NPGs,genesList$hk,genesList$PNGs),c(genesList$NPGs,genesList$hk,genesList$PNGs)]

f1 = colorRamp2(seq(-0.5,0.5, length = 3), c("#DC0000B2", "white","#3C5488B2" ))

split.genes <- factor(c(rep("NPGs",length(genesList[["NPGs"]])),
                         rep("HK",length(genesList[["hk"]])),
                         rep("PNGs",length(genesList[["PNGs"]]))
                        ),
                         levels = c("NPGs","HK","PNGs"))

lgd = Legend(col_fun = f1, title = "COTAN coex")

htmp <- Heatmap(as.matrix(coexMat),
        #width = ncol(coexMat)*unit(2.5, "mm"), 
        height = nrow(coexMat)*unit(3, "mm"),
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        col = f1,
        row_names_side = "left",
        row_names_gp = gpar(fontsize = 11),
        column_names_gp  = gpar(fontsize = 11),
        column_split = split.genes,
        row_split = split.genes,
        cluster_row_slices = FALSE, 
    cluster_column_slices = FALSE,
    heatmap_legend_param = list(
        title = "COTAN coex", at = c(-0.5, 0, 0.5),direction = "horizontal",
        labels = c("-0.5", "0", "0.5")
    )
   )
draw(htmp, heatmap_legend_side="right")
```

```{r, echo=FALSE}
pvalueMat.big <- calculatePValue(obj)#, geneSubsetCol = int.genes, geneSubsetRow = int.genes)
```

```{r, echo=FALSE}

write_fst(as.data.frame(coexMat.big), paste0("CoexData/CotanCoex_",file_code,
                              ".fst"),compress = 100)
write_fst(as.data.frame(pvalueMat.big), paste0("CoexData/CotanPValues_",
                              file_code,".fst"),compress = 100)


rm(coexMat.big)
rm(pvalueMat.big)
```

```{r}
GDI_DF <- calculateGDI(obj)
GDI_DF$geneType <- NA
for (cat in names(genesList)) {
  GDI_DF[rownames(GDI_DF) %in% genesList[[cat]],]$geneType <- cat
}

GDI_DF$GDI_centered <- scale(GDI_DF$GDI,center = T,scale = T)

write.csv(GDI_DF,paste0("CoexData/","Variance_GDI_genes",file_code,".csv"))


GDI_DF[unlist(genesList),]
```

```{r}
GDIPlot(obj,GDIIn = GDI_DF, genes = genesList,GDIThreshold = 1.4)
```

## Seurat correlation

```{r}
srat<- CreateSeuratObject(counts = getRawData(obj), 
                          project = project, 
                          min.cells = 3, 
                          min.features = 200)
srat[["percent.mt"]] <- PercentageFeatureSet(srat, pattern = "^mt-")
srat <- NormalizeData(srat)
srat <- FindVariableFeatures(srat, selection.method = "vst", nfeatures = 2000)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(srat)

plot1$data$centered_variance <- scale(plot1$data$variance.standardized,
                                      center = T,scale = F)

write.csv(plot1$data,paste0("CoexData/",
                            "Variance_Seurat_genes",
                            getMetadataElement(obj, 
                                               datasetTags()[["cond"]]),".csv"))

LabelPoints(plot = plot1, points = c(genesList$NPGs,genesList$PNGs,genesList$layers), repel = TRUE)

```

```{r}
LabelPoints(plot = plot1, points = c(genesList$hk), repel = TRUE)
```

```{r}
all.genes <- rownames(srat)
srat <- ScaleData(srat, features = all.genes)
seurat.data = GetAssayData(srat[["RNA"]],layer = "data")
```

```{r}
corr.pval.list <- correlation_pvalues(data = seurat.data,
                                      int.genes,
                                      n.cells = getNumCells(obj))

seurat.data.cor.big <- as.matrix(Matrix::forceSymmetric(corr.pval.list$data.cor, uplo = "U"))

htmp <- correlation_plot(seurat.data.cor.big, 
                         genesList, title="Seurat corr")


p_values.fromSeurat <- corr.pval.list$p_values
seurat.data.cor.big <- corr.pval.list$data.cor

rm(corr.pval.list)
gc()

draw(htmp, heatmap_legend_side="right")
```

```{r, echo=FALSE}

write_fst(as.data.frame(p_values.fromSeurat),path = paste0("CoexData/SeuratPValues_",file_code,".fst"),compress = 100)

write_fst(as.data.frame(seurat.data.cor.big),
          path = paste0("CoexData/SeuratCorr_",
                              file_code,
                              ".fst"), compress = 100)

rm(seurat.data.cor.big)
rm(p_values.fromSeurat)
```

### Seurat SC Transform

```{r}
srat <-  SCTransform(srat, 
                     method = "glmGamPoi", 
                     vars.to.regress = "percent.mt", 
                     verbose = FALSE)

seurat.data <- GetAssayData(srat[["SCT"]],layer = "data")

#Remove genes with all zeros
seurat.data <-seurat.data[rowSums(seurat.data) > 0,]


corr.pval.list <- correlation_pvalues(seurat.data,
                                      int.genes,
                                      n.cells = getNumCells(obj))

seurat.data.cor.big <- as.matrix(Matrix::forceSymmetric(corr.pval.list$data.cor, uplo = "U"))

htmp <- correlation_plot(seurat.data.cor.big, 
                         genesList, title="Seurat corr SCT")



p_values.fromSeurat <- corr.pval.list$p_values
seurat.data.cor.big <- corr.pval.list$data.cor

rm(corr.pval.list)
gc()


draw(htmp, heatmap_legend_side="right")
```

```{r,echo=FALSE}
plot1 <- VariableFeaturePlot(srat)

plot1$data$centered_variance <- scale(plot1$data$residual_variance,
                                      center = T,scale = F)

write.csv(plot1$data,paste0("CoexData/",
                            "Variance_SeuratSCT_genes",
                            getMetadataElement(obj, 
                                               datasetTags()[["cond"]]),".csv"))

write_fst(as.data.frame(seurat.data.cor.big),path = paste0("CoexData/SeuratCorrSCT_",file_code,".fst"), compress = 100)

write_fst(as.data.frame(p_values.fromSeurat),path =  paste0("CoexData/SeuratPValuesSCT_", file_code,".fst"))


rm(seurat.data.cor.big)
rm(p_values.fromSeurat)
```

## Monocle

```{r}
library(monocle3)
```

```{r}
cds <- new_cell_data_set(getRawData(obj),
                         cell_metadata = getMetadataCells(obj),
                         gene_metadata = getMetadataGenes(obj)
                         )
cds <- preprocess_cds(cds, num_dim = 100)

normalized_counts <- normalized_counts(cds)
```

```{r}
#Remove genes with all zeros
normalized_counts <- normalized_counts[rowSums(normalized_counts) > 0,]


corr.pval.list <- correlation_pvalues(normalized_counts,
                                      int.genes,
                                      n.cells = getNumCells(obj))

rm(normalized_counts)

monocle.data.cor.big <- as.matrix(Matrix::forceSymmetric(corr.pval.list$data.cor, uplo = "U"))

htmp <- correlation_plot(data.cor.big = monocle.data.cor.big,
                         genesList,
                         title = "Monocle corr")


p_values.from.monocle <- corr.pval.list$p_values
monocle.data.cor.big <- corr.pval.list$data.cor

rm(corr.pval.list)
gc()

draw(htmp, heatmap_legend_side="right")
```

```{r,echo=FALSE}
write_fst(as.data.frame(monocle.data.cor.big),
       path = paste0("CoexData/MonocleCorr_",
                              file_code,
                              ".fst"),compress = 100)
write_fst(as.data.frame(p_values.from.monocle),
       path = paste0("CoexData/MonoclePValues_",
                              file_code,".fst"),compress = 100)



rm(p_values.from.monocle)
rm(monocle.data.cor.big)
```

## ScanPy

```{r}
library(reticulate)

dirOutScP <- paste0("CoexData/ScanPy/")
if (!dir.exists(dirOutScP)) {
  dir.create(dirOutScP)
}

Sys.setenv(RETICULATE_PYTHON = "../../../bin/python3")
  py <- import("sys")
  
  source_python("src/scanpyGenesExpression.py")
  scanpyFDR(getRawData(obj), 
                  getMetadataCells(obj), 
                  getMetadataGenes(obj), 
                   "mt", 
                   dirOutScP, 
                   file_code,
            int.genes)

normalized_counts <- read.csv(paste0(dirOutScP,
                                     file_code,"_Scampy_expression_all_genes.gz"),header = T,row.names = 1)

normalized_counts <- t(normalized_counts)

```

```{r}
#Remove genes with all zeros
normalized_counts <-normalized_counts[rowSums(normalized_counts) > 0,]


corr.pval.list <- correlation_pvalues(normalized_counts,
                                      int.genes,
                                      n.cells = getNumCells(obj))

ScanPy.data.cor.big <- as.matrix(Matrix::forceSymmetric(corr.pval.list$data.cor, uplo = "U"))


htmp <- correlation_plot(data.cor.big = ScanPy.data.cor.big,
                         genesList,
                         title = "ScanPy corr")

p_values.from.ScanPy <- corr.pval.list$p_values
ScanPy.data.cor.big <- corr.pval.list$data.cor

rm(corr.pval.list)
gc()

draw(htmp, heatmap_legend_side="right")
```

```{r,echo=FALSE}
write_fst(as.data.frame(ScanPy.data.cor.big), path = paste0("CoexData/ScanPyCorr_", file_code,      ".fst"),compress = 100)
write_fst(as.data.frame(p_values.from.ScanPy),
       path = paste0("CoexData/ScanPyPValues_",
                              file_code,".fst"), compress = 100)

rm(p_values.fromScanPy)
rm(ScanPy.data.cor.big)
```

------------------------------------------------------------------------

```{r}
Sys.time()
```

```{r}
sessionInfo()
```