Spearman Correlation

Setup knitr and load utility functions

knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir="E:/DISC/reproducibility")

utilities_path = "./source/utilities.r"
source(utilities_path)

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
reldist: Relative Distribution Methods
Version 1.6-6 created on 2016-10-07.
copyright (c) 2003, Mark S. Handcock, University of California-Los Angeles
 For citation information, type citation("reldist").
 Type help(package="reldist") to get started.

Settings

method_name = c("Raw", "DISC", "scVI", "MAGIC", "DCA", "scScope", "DeepImpute", "VIPER", "scImpute")
method_color = c("#A5A5A5", "#E83828", "#278BC4", "#EADE36", "#198B41", "#920783", "#F8B62D", "#8E5E32", "#1E2B68")
text_color = rep("black", length(method_name))
names(text_color) = method_name
text_color["DISC"] = "red"

Here, we use Spearman Correlation Coefficient as an indicator to evaluate the similarity between bulk data and imputed single data. For replicated bulk samples, we average their profiles as this paper (https://www.biorxiv.org/content/10.1101/2020.01.29.925974v1).

gene_bulk_fpkm = readh5_loom("./data/JURKAT_293T/bulk_fpkm.loom")

[1] "./data/JURKAT_293T/bulk_fpkm.loom"
[1] 58011     4

gene_bulk_fpkm_log1p = log1p(gene_bulk_fpkm)
bulk_data = cbind(rowMeans(gene_bulk_fpkm_log1p[, 1:2]), rowMeans(gene_bulk_fpkm_log1p[, 3:4]))
colnames(bulk_data) = c("293T", "Jurkat")

cor_list = list()

293T

gene_bc_mat = readh5_loom("./data/JURKAT_293T/raw_293T.loom")

[1] "./data/JURKAT_293T/raw_293T.loom"
[1] 32738  2885

gene_filter = gene_selection(gene_bc_mat, 10)
sc_gene = rownames(gene_bc_mat)[gene_filter]
compare_gene = intersect(sc_gene, rownames(bulk_data))
cor_list[["293T"]] = c()
for(ii in method_name){
  switch(ii, 
         "Raw" = {
           gene_bc_filt = gene_bc_mat[sc_gene, ]
         },
         "DISC" = {
           gene_bc_mat = readh5_loom(paste0("./data/JURKAT_293T/DISC_293T.loom"))
           gene_bc_filt = gene_bc_mat[sc_gene, ]
         }, {
           gene_bc_mat = readh5_imputation(paste0("./data/JURKAT_293T/", ii, "_293T.hdf5"))
           gene_bc_filt = gene_bc_mat[sc_gene, ]
         }
  )
  seurat_obj = CreateSeuratObject(counts = as.data.frame(gene_bc_filt), min.cells = 0, min.features = 0)
  seurat_obj = NormalizeData(seurat_obj, normalization.method = "LogNormalize", scale.factor = 10000)
  sc_norm = as.matrix(seurat_obj[["RNA"]]@data)
  rownames(sc_norm) = sc_gene
  sc_data = sc_norm[compare_gene, ]
  cor_list[["293T"]] = c(cor_list[["293T"]], cor(rowMeans(sc_data), bulk_data[compare_gene, 1], method = "spearman"))
}

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/DISC_293T.loom"
[1] 32738  2885

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/scVI_293T.hdf5"
[1] 11974  2885

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/MAGIC_293T.hdf5"
[1] 11974  2885

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/DCA_293T.hdf5"
[1] 11974  2885

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/scScope_293T.hdf5"
[1] 11974  2885

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/DeepImpute_293T.hdf5"
[1] 11974  2885

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/VIPER_293T.hdf5"
[1] 11974  2885

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/scImpute_293T.hdf5"
[1] 11974  2885

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

names(cor_list[["293T"]]) = method_name

Jurkat

gene_bc_mat = readh5_loom("./data/JURKAT_293T/raw_Jurkat.loom")

[1] "./data/JURKAT_293T/raw_Jurkat.loom"
[1] 32738  3258

gene_filter = gene_selection(gene_bc_mat, 10)
sc_gene = rownames(gene_bc_mat)[gene_filter]
compare_gene = intersect(sc_gene, rownames(bulk_data))
cor_list[["Jurkat"]] = c()
for(ii in method_name){
  switch(ii, 
         "Raw" = {
           gene_bc_filt = gene_bc_mat[sc_gene, ]
         },
         "DISC" = {
           gene_bc_mat = readh5_loom(paste0("./data/JURKAT_293T/DISC_Jurkat.loom"))
           gene_bc_filt = gene_bc_mat[sc_gene, ]
         }, {
           gene_bc_mat = readh5_imputation(paste0("./data/JURKAT_293T/", ii, "_Jurkat.hdf5"))
           gene_bc_filt = gene_bc_mat[sc_gene, ]
         }
  )
  seurat_obj = CreateSeuratObject(counts = as.data.frame(gene_bc_filt), min.cells = 0, min.features = 0)
  seurat_obj = NormalizeData(seurat_obj, normalization.method = "LogNormalize", scale.factor = 10000)
  sc_norm = as.matrix(seurat_obj[["RNA"]]@data)
  rownames(sc_norm) = sc_gene
  sc_data = sc_norm[compare_gene, ]
  cor_list[["Jurkat"]] = c(cor_list[["Jurkat"]], cor(rowMeans(sc_data), bulk_data[compare_gene, 2], method = "spearman"))
}

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/DISC_Jurkat.loom"
[1] 32738  3258

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/scVI_Jurkat.hdf5"
[1] 11293  3258

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/MAGIC_Jurkat.hdf5"
[1] 11293  3258

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/DCA_Jurkat.hdf5"
[1] 11293  3258

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/scScope_Jurkat.hdf5"
[1] 11293  3258

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/DeepImpute_Jurkat.hdf5"
[1] 11293  3258

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/VIPER_Jurkat.hdf5"
[1] 11293  3258

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/scImpute_Jurkat.hdf5"
[1] 11293  3258

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

names(cor_list[["Jurkat"]]) = method_name

JURKAT_293T

gene_bc_mat = readh5_loom("./data/JURKAT_293T/raw.loom")

[1] "./data/JURKAT_293T/raw.loom"
[1] 32738  6143

gene_filter = gene_selection(gene_bc_mat, 10)
sc_gene = rownames(gene_bc_mat)[gene_filter]
cell_Jurkat = colnames(gene_bc_mat)[grep("JURKAT", colnames(gene_bc_mat))]
cell_293T = colnames(gene_bc_mat)[grep("293T", colnames(gene_bc_mat))]
compare_gene = intersect(sc_gene, rownames(bulk_data))
cor_list[["JURKAT_293T"]] = c()
for(ii in method_name){
  if(ii == "Raw"){
    gene_bc_filt = gene_bc_mat[sc_gene, ]
  }else{
    gene_bc_mat = readh5_loom(paste0("./data/JURKAT_293T/", ii, "_JURKAT_293T.loom"))
    gene_bc_filt = gene_bc_mat[sc_gene, ]
  }
  seurat_obj = CreateSeuratObject(counts = as.data.frame(gene_bc_filt), min.cells = 0, min.features = 0)
  seurat_obj = NormalizeData(seurat_obj, normalization.method = "LogNormalize", scale.factor = 10000)
  sc_norm = as.matrix(seurat_obj[["RNA"]]@data)
  rownames(sc_norm) = sc_gene
  sc_data = sc_norm[compare_gene, ]
  cor_list[["JURKAT_293T"]] = c(cor_list[["JURKAT_293T"]], cor(rowMeans(sc_data[, cell_293T]) - rowMeans(sc_data[, cell_Jurkat]), bulk_data[compare_gene, 1] - bulk_data[compare_gene, 2], method = "spearman"))
}

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/DISC_JURKAT_293T.loom"
[1] 32738  6143

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/scVI_JURKAT_293T.loom"
[1] 13328  6143

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/MAGIC_JURKAT_293T.loom"
[1] 13328  6143

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/DCA_JURKAT_293T.loom"
[1] 13328  6143

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/scScope_JURKAT_293T.loom"
[1] 13328  6143

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/DeepImpute_JURKAT_293T.loom"
[1] 13328  6143

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/VIPER_JURKAT_293T.loom"
[1] 13328  6143

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "./data/JURKAT_293T/scImpute_JURKAT_293T.loom"
[1] 13328  6143

Feature names cannot have underscores ('_'), replacing with dashes ('-')Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

names(cor_list[["JURKAT_293T"]]) = method_name

Results

barplot_usage(cor_list[["Jurkat"]], main = "", cex.main = 1.5, bar_color = method_color, text_color = text_color, use_data_order = T, ylab = "Jurkat", cex.lab = 1.5, font.main = 1, ylim = c(-0.1, 1), use_border = F, decreasing = T)

barplot_usage(cor_list[["293T"]], main = "", cex.main = 1.5, bar_color = method_color, text_color = text_color, use_data_order = T, ylab = "293T", cex.lab = 1.5, font.main = 1, ylim = c(-0.1, 1), use_border = F, decreasing = T)

barplot_usage(cor_list[["JURKAT_293T"]], main = "", cex.main = 1.5, bar_color = method_color, text_color = text_color, use_data_order = T, ylab = "JURKAT_293T", cex.lab = 1.5, font.main = 1, ylim = c(-0.1, 1), use_border = F, decreasing = T)