Seurat correlation vs cotan heatmap

# Comparition between COTAN and Seurat
library(dplyr)
#source("cotan_output_functions.R")
library(Seurat)
library(Matrix)
library(fitdistrplus)

Loading required package: survival

library("ggsci")
library("ggplot2")
library(ggrepel)
library("gridExtra")

Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine

library(ggpubr)

Attaching package: ‘ggpubr’

The following object is masked from ‘package:dendextend’:

    rotate

library(COTAN)

hk = c("Calm1","Cox6b1","Ppia","Rpl18","Cox7c","Erh","H3f3a","Taf1","Taf2","Gapdh","Actb","Golph3",
       "Mtmr12","Zfr","Sub1","Tars","Amacr")
tf1 = c("Nes","Vim","Sox2","Sox1","Notch1", "Hes1","Hes5","Pax6")
tf2 = c("Map2","Tubb3","Neurod1","Nefm","Nefl","Dcx","Tbr1")
layer = c("Reln","Lhx5","Cux1","Satb2","Tle1","Mef2c","Rorb","Sox5","Bcl11b","Fezf2","Foxp2")

input_dir = "Data/mouse_dentate_gyrus/"
 
E16.5_cotan = readRDS("Data/E16.5_cleaned.cotan.RDS")

t = "E16.5"

pbmc.data  =  read.csv("Data/mouse_dentate_gyrus/separation_age/E16.5.csv", header = T, row.names = 1)

#cells = read.csv(paste(input_dir,"cells_",t,".csv", sep = ""), header = T, row.names = 1)

#out_dir = "../results/2019.12.14/"
#load_files(imput_dir = input_dir,t="E16.5",s="H1")
#colnames(pbmc.data)=paste("E",t,"_",colnames(pbmc.data), sep = "")
pbmc.data = pbmc.data[, colnames(E16.5_cotan@raw)]

rownames(pbmc.data)[grep("_",rownames(pbmc.data))] =  gsub("_","-",rownames(pbmc.data)[grep("_",rownames(pbmc.data))]) 
#colnames(pbmc.data)=stringr::str_remove(colnames(pbmc.data), "\\.$")
#colnames(pbmc.data)[grep(".",colnames(pbmc.data))] =  gsub("."," ",colnames(pbmc.data)[grep(".",colnames(pbmc.data))]) 
pbmc.data = as(as.matrix(pbmc.data), "sparseMatrix")
lib.size = colSums(pbmc.data)
hist(lib.size, prob=TRUE, breaks=145)

fit <- fitdist(lib.size, "nbinom")

NaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs produced

prob.treshold = qnbinom((1-0.001), size=fit$estimate[1], mu=fit$estimate[2])
pbmc.data = pbmc.data[,names(lib.size[lib.size<prob.treshold])]
lib.size = colSums(pbmc.data)
hist(lib.size, prob=TRUE, breaks=145)

pbmc <- CreateSeuratObject(counts = pbmc.data, project = "hipp_1_E16.5", min.cells = 3, min.features = 200)
pbmc

An object of class Seurat 
15508 features across 2228 samples within 1 assay 
Active assay: RNA (15508 features, 0 variable features)

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^mt-")
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))

CombinePlots is being deprecated. Plots should now be combined using the patchwork system.

subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 4000 & percent.mt < 7.5)

An object of class Seurat 
15508 features across 2080 samples within 1 assay 
Active assay: RNA (15508 features, 0 variable features)

length(pbmc$nFeature_RNA)

[1] 2228

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))

CombinePlots is being deprecated. Plots should now be combined using the patchwork system.

pbmc <- NormalizeData(pbmc)

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

pbmc[["RNA"]]@data[1:10,1:10]

10 x 10 sparse Matrix of class "dgCMatrix"

   [[ suppressing 10 column names ‘X10X82_1_CAGAACAAGTTCTG.’, ‘X10X82_1_TCACTTCTCGCATC.’, ‘X10X82_1_TCACTCATATGCTG.’ ... ]]

                                                                                                            
0610007P14Rik 1.7126992 1.143042 1.077085 0.9062106 2.082862 0.8906782 1.3936138 .        1.680767 1.8883605
0610009B22Rik 0.7589239 1.143042 .        0.9062106 1.206710 0.8906782 .         1.700242 .        .        
0610009L18Rik .         .        .        0.9062106 .        .         .         .        .        .        
0610009O20Rik 1.1853870 .        .        .         .        .         .         .        .        .        
0610010F05Rik .         1.143042 .        .         .        0.8906782 0.9221505 .        .        .        
0610010K14Rik 1.7126992 .        .        0.9062106 .        .         .         .        .        0.8763564
0610011F06Rik .         .        1.077085 0.9062106 1.206710 .         0.9221505 .        .        0.8763564
0610012G03Rik 1.7126992 1.143042 1.918146 1.3736782 1.206710 1.3541748 0.9221505 2.042896 1.158401 1.3361221
0610025J13Rik .         .        .        .         .        .         0.9221505 .        .        .        
0610030E20Rik .         .        .        .         .        0.8906782 .         .        .        .        

Code to produce the heatmap of Seurat correlation matrix vs cotan heatmap on selected genes.

seurat.data = pbmc[["RNA"]]@data
seurat.data.cor.pearson = cor(t(as.matrix(seurat.data)), method = "pearson")
#write.csv(seurat.data.cor, file = "../results/2019.12.16/E16.5_hipp_Seurat_correlations.csv" )

partial.coex.pearson = seurat.data.cor.pearson[rownames(seurat.data.cor.pearson) %in% c(tf1,tf2,hk),colnames(seurat.data.cor.pearson) %in% c(tf1,tf2,hk)]
diag(partial.coex.pearson) = 0

seurat.data.cor.spearman = cor(t(as.matrix(seurat.data)), method = "spearman")
#write.csv(seurat.data.cor, file = "../results/2019.12.16/E16.5_hipp_Seurat_correlations.csv" )

partial.coex.spearman = seurat.data.cor.spearman[rownames(seurat.data.cor.spearman) %in% c(tf1,tf2,hk),colnames(seurat.data.cor.spearman) %in% c(tf1,tf2,hk)]
diag(partial.coex.spearman) = 0

1 Pearson correlation

partial.coex.pearson = reshape2::melt(partial.coex.pearson)
colnames(partial.coex.pearson) = c("g1","g2","corr")

partial.coex.pearson$g1 <- factor(partial.coex.pearson$g1, c(tf1,hk,tf2))
partial.coex.pearson$g2 <- factor(partial.coex.pearson$g2, c(tf1,hk,tf2))

  P = ggplot(partial.coex.pearson) + 
  geom_tile(aes(x=g1,y=g2, fill = corr),colour = "black", show.legend = TRUE) +
#  facet_grid( g1 ~ g2  ,scales = "free", space = "free") + 
    scale_fill_gradient2(mid = "white",limits=c(-1, 1),low = "#DC0000B2", high = "#3C5488B2")+
  #scale_fill_gradient2(low = "darkred", mid = "white",  high = "darkblue", midpoint = 0,na.value = "grey80", space = "Lab", guide = "colourbar", aesthetics = "fill", limits = lim_coex, oob=scales::squish)+ theme(legend.position="bottom")+
  theme(#legend.title = element_blank(),
    #strip.text.x = element_text(color = "red"),
    #axis.text.y = element_text(color = ),
    axis.text.x = element_text(angle=45,hjust=1,vjust=1.0),
    legend.position="bottom"
  ) #+geom_text(aes(label=ifelse(t_hk == "hk", "H","")), color="grey", size=3)
  P

2 Spearman correlation

partial.coex.spearman = reshape2::melt(partial.coex.spearman)
colnames(partial.coex.spearman) = c("g1","g2","corr")

partial.coex.spearman$g1 <- factor(partial.coex.spearman$g1, c(tf1,hk,tf2))
partial.coex.spearman$g2 <- factor(partial.coex.spearman$g2, c(tf1,hk,tf2))

  S = ggplot(partial.coex.spearman) + 
  geom_tile(aes(x=g1,y=g2, fill = corr),colour = "black", show.legend = TRUE) +
#  facet_grid( g1 ~ g2  ,scales = "free", space = "free") + 
    scale_fill_gradient2(mid = "white",limits=c(-1, 1),low = "#DC0000B2", high = "#3C5488B2")+
  #scale_fill_gradient2(low = "darkred", mid = "white",  high = "darkblue", midpoint = 0,na.value = "grey80", space = "Lab", guide = "colourbar", aesthetics = "fill", limits = lim_coex, oob=scales::squish)+ theme(legend.position="bottom")+
  theme(#legend.title = element_blank(),
    #strip.text.x = element_text(color = "red"),
    #axis.text.y = element_text(color = ),
    axis.text.x = element_text(angle=45,hjust=1,vjust=1.0),
    legend.position="bottom"
  ) #+geom_text(aes(label=ifelse(t_hk == "hk", "H","")), color="grey", size=3)
  S

3 COTAN coex

#load_data3.0(input_dir, cond = t, genes = c(tf1,tf2,hk), prefix = "p_value_")

# COTAN: after loading in memeory the coex matrix 
p_value_E16.5 = get.pval(object = E16.5_cotan, gene.set.col = c(tf1,tf2,hk),
                         gene.set.row = c(tf1,tf2,hk))

 [1] "Nes"     "Vim"     "Sox2"    "Sox1"    "Notch1"  "Hes1"    "Hes5"    "Pax6"    "Map2"    "Tubb3"  
[11] "Neurod1" "Nefm"    "Nefl"    "Dcx"     "Tbr1"    "Calm1"   "Cox6b1"  "Ppia"    "Rpl18"   "Cox7c"  
[21] "Erh"     "H3f3a"   "Taf1"    "Taf2"    "Gapdh"   "Actb"    "Golph3"  "Mtmr12"  "Zfr"     "Sub1"   
[31] "Tars"    "Amacr"  
[1] "Get p-values on a set of genes on columns on a set of genes on rows"
[1] "Using function S"
[1] "function to generate S "

partial.coex.cotan = E16.5_cotan@coex[rownames(E16.5_cotan@coex) %in% c(tf1,tf2,hk),colnames(E16.5_cotan@coex) %in% c(tf1,tf2,hk)]

#partial.pval.cotan = p_value_E16.5[rownames(p_value_E16.5) %in% c(tf1,tf2,hk),colnames(p_value_E16.5) %in% c(tf1,tf2,hk)]
partial.pval.cotan =p_value_E16.5
#partial.pval.cotan = partial.pval.cotan <= 0.05
#partial.coex.cotan[!partial.pval.cotan] <- 0

  partial.coex.cotan = reshape2::melt(as.matrix(partial.coex.cotan))
  colnames(partial.coex.cotan) = c("g1","g2","coex")
  for (n in c(1:nrow(partial.coex.cotan))) {
    if (partial.coex.cotan[n,"g1"] == partial.coex.cotan[n,"g2"]) {
      partial.coex.cotan[n,"coex"]=0
    }
      
  }
  
  partial.coex.cotan$g1 <- factor(partial.coex.cotan$g1, c(tf1,hk,tf2))
  partial.coex.cotan$g2 <- factor(partial.coex.cotan$g2, c(tf1,hk,tf2))
  
C = ggplot(partial.coex.cotan) + 
    geom_tile(aes(x=g1,y=g2, fill = coex),colour = "black", show.legend = TRUE) +
    #  facet_grid( g1 ~ g2  ,scales = "free", space = "free") + 
    scale_fill_gradient2(mid = "white",limits=c(round(min(partial.coex.cotan$coex),digits = 0), round(max(partial.coex.cotan$coex),digits = 0)),low = "#DC0000B2", high = "#3C5488B2")+
    #scale_fill_gradient2(low = "darkred", mid = "white",  high = "darkblue", midpoint = 0,na.value = "grey80", space = "Lab", guide = "colourbar", aesthetics = "fill", limits = lim_coex, oob=scales::squish)+ theme(legend.position="bottom")+
    theme(#legend.title = element_blank(),
      #strip.text.x = element_text(color = "red"),
      #axis.text.y = element_text(color = ),
      axis.text.x = element_text(angle=45,hjust=1,vjust=1.0),
      legend.position="bottom"
    )#+geom_text(aes(label=ifelse(t_hk == "hk", "H","")), color="grey", size=3)
  

#figure <- ggarrange(C, S,
 #                   labels = c("Co.", "Sp."),
 #                   ncol = 2, nrow = 1)
#figure
C

figure <- ggarrange(C, P,
                    labels = c("Co.", "Pe."),
                    ncol = 2, nrow = 1)
figure

sessionInfo()

R version 4.0.4 (2021-02-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggpubr_0.4.0          gridExtra_2.3         ggsci_2.9             fitdistrplus_1.1-3   
 [5] survival_3.2-7        gplots_3.1.1          patchwork_1.1.1       SeuratObject_4.0.0   
 [9] Seurat_4.0.0          cluster_2.1.1         WGCNA_1.70-3          fastcluster_1.1.25   
[13] dynamicTreeCut_1.63-1 dendextend_1.14.0     MASS_7.3-53.1         htmlwidgets_1.5.3    
[17] forcats_0.5.1         stringr_1.4.0         dplyr_1.0.4           purrr_0.3.4          
[21] readr_1.4.0           tidyr_1.1.2           tibble_3.0.6          tidyverse_1.3.0      
[25] plotly_4.9.3          Rtsne_0.15            factoextra_1.0.7      ggrepel_0.9.1        
[29] ggplot2_3.3.3         Matrix_1.3-2          data.table_1.13.6     COTAN_0.1.0          

loaded via a namespace (and not attached):
  [1] reticulate_1.18       R.utils_2.10.1        tidyselect_1.1.0      RSQLite_2.2.3        
  [5] AnnotationDbi_1.52.0  grid_4.0.4            munsell_0.5.0         preprocessCore_1.52.1
  [9] codetools_0.2-18      ica_1.0-2             future_1.21.0         miniUI_0.1.1.1       
 [13] withr_2.4.1           colorspace_2.0-0      Biobase_2.50.0        filelock_1.0.2       
 [17] knitr_1.31            rstudioapi_0.13       stats4_4.0.4          ROCR_1.0-11          
 [21] ggsignif_0.6.1        tensor_1.5            listenv_0.8.0         labeling_0.4.2       
 [25] polyclip_1.10-0       bit64_4.0.5           farver_2.0.3          basilisk_1.2.1       
 [29] parallelly_1.23.0     vctrs_0.3.6           generics_0.1.0        xfun_0.20            
 [33] R6_2.5.0              doParallel_1.0.16     clue_0.3-58           bitops_1.0-6         
 [37] cachem_1.0.3          spatstat.utils_2.0-0  assertthat_0.2.1      promises_1.1.1       
 [41] scales_1.1.1          nnet_7.3-15           gtable_0.3.0          Cairo_1.5-12.2       
 [45] globals_0.14.0        goftest_1.2-2         rlang_0.4.10          GlobalOptions_0.1.2  
 [49] splines_4.0.4         rstatix_0.7.0         lazyeval_0.2.2        impute_1.64.0        
 [53] checkmate_2.0.0       broom_0.7.5           yaml_2.2.1            reshape2_1.4.4       
 [57] abind_1.4-5           modelr_0.1.8          crosstalk_1.1.1       backports_1.2.1      
 [61] httpuv_1.5.5          Hmisc_4.5-0           tools_4.0.4           ellipsis_0.3.1       
 [65] jquerylib_0.1.3       RColorBrewer_1.1-2    BiocGenerics_0.36.0   ggridges_0.5.3       
 [69] latex2exp_0.4.0       Rcpp_1.0.6            plyr_1.8.6            base64enc_0.1-3      
 [73] basilisk.utils_1.2.2  rpart_4.1-15          deldir_0.2-10         pbapply_1.4-3        
 [77] GetoptLong_1.0.5      viridis_0.5.1         cowplot_1.1.1         S4Vectors_0.28.1     
 [81] zoo_1.8-8             haven_2.3.1           fs_1.5.0              magrittr_2.0.1       
 [85] scattermore_0.7       openxlsx_4.2.3        circlize_0.4.12       lmtest_0.9-38        
 [89] reprex_1.0.0          RANN_2.6.1            matrixStats_0.58.0    hms_1.0.0            
 [93] mime_0.9              evaluate_0.14         xtable_1.8-4          jpeg_0.1-8.1         
 [97] rio_0.5.16            readxl_1.3.1          IRanges_2.24.1        shape_1.4.5          
[101] compiler_4.0.4        KernSmooth_2.23-18    crayon_1.4.0          R.oo_1.24.0          
[105] htmltools_0.5.1.1     mgcv_1.8-33           later_1.1.0.1         Formula_1.2-4        
[109] lubridate_1.7.9.2     DBI_1.1.1             dbplyr_2.1.0          ComplexHeatmap_2.6.2 
[113] rappdirs_0.3.3        car_3.0-10            cli_2.3.0             R.methodsS3_1.8.1    
[117] parallel_4.0.4        igraph_1.2.6          pkgconfig_2.0.3       foreign_0.8-81       
[121] xml2_1.3.2            foreach_1.5.1         bslib_0.2.4           rvest_0.3.6          
[125] digest_0.6.27         sctransform_0.3.2     RcppAnnoy_0.0.18      spatstat.data_2.0-0  
[129] rmarkdown_2.7         cellranger_1.1.0      leiden_0.3.7          htmlTable_2.1.0      
[133] uwot_0.1.10           curl_4.3              gtools_3.8.2          shiny_1.6.0          
[137] rjson_0.2.20          lifecycle_0.2.0       nlme_3.1-152          jsonlite_1.7.2       
[141] carData_3.0-4         viridisLite_0.3.0     pillar_1.4.7          lattice_0.20-41      
[145] fastmap_1.1.0         httr_1.4.2            GO.db_3.12.1          glue_1.4.2           
[149] zip_2.1.1             spatstat_1.64-1       png_0.1-7             iterators_1.0.13     
[153] bit_4.0.4             stringi_1.5.3         sass_0.3.1            blob_1.2.1           
[157] caTools_1.18.1        latticeExtra_0.6-29   memoise_2.0.0         irlba_2.3.3          
[161] future.apply_1.7.0