Dropout event recovery

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.

Load raw data and downsampling

Here, we use MELANOMA dataset as an example.
The fraction of sampling is set to 50% of the original library size across cells.
Only 1 replicate will be generated here as an example.
Only imputed genes will be kept for comparison.

raw_data = readh5_loom("./data/MELANOMA/raw.loom")

[1] "./data/MELANOMA/raw.loom"
[1] 32287  8498

ds_dir = "./data/MELANOMA/ds_0.5/r1"
dir.create(ds_dir, showWarnings = F, recursive = T)
output_dir = "./results/MELANOMA/ds_0.5/r1"
dir.create(output_dir, showWarnings = F, recursive = T)
observed_path = paste0(ds_dir, "/observed.loom")
if(!file.exists(observed_path)){
  observed_data = downsampling_cell(0.5, raw_data)
  save_h5(observed_path, t(observed_data))
}else{
  observed_data = readh5_loom(observed_path)
}

[1] "./data/MELANOMA/ds_0.5/r1/observed.loom"
[1] 32287  8498

gene_filter = gene_selection(observed_data, 10)
raw_data = raw_data[gene_filter, ]
observed_data = observed_data[gene_filter, ]
used_genes = rownames(observed_data)
print(dim(raw_data))

[1] 12509  8498

print(dim(observed_data))

[1] 12509  8498

data_list = list(Raw = raw_data, Observed = observed_data)
rm(raw_data, observed_data)

Load downsampling data and imputation results

We run imputation using downsampling data.
Here, we load all imputation results.

data_list[["DISC"]] = readh5_loom(paste0(ds_dir, "/DISC.loom"))[used_genes, ]

[1] "./data/MELANOMA/ds_0.5/r1/DISC.loom"
[1] 32287  8498

print(dim(data_list[["DISC"]]))

[1] 12509  8498

data_list[["scVI"]] = readh5_imputation(paste0(ds_dir, "/scVI.hdf5"))

[1] "./data/MELANOMA/ds_0.5/r1/scVI.hdf5"
[1] 12509  8498

data_list[["MAGIC"]] = readh5_imputation(paste0(ds_dir, "/MAGIC.hdf5"))

[1] "./data/MELANOMA/ds_0.5/r1/MAGIC.hdf5"
[1] 12509  8498

data_list[["DCA"]] = readh5_imputation(paste0(ds_dir, "/DCA.hdf5"))

[1] "./data/MELANOMA/ds_0.5/r1/DCA.hdf5"
[1] 12509  8498

data_list[["scScope"]] = readh5_imputation(paste0(ds_dir, "/scScope.hdf5"))

[1] "./data/MELANOMA/ds_0.5/r1/scScope.hdf5"
[1] 12509  8498

data_list[["DeepImpute"]] = readh5_imputation(paste0(ds_dir, "/DeepImpute.hdf5"))

[1] "./data/MELANOMA/ds_0.5/r1/DeepImpute.hdf5"
[1] 12509  8498

data_list[["VIPER"]] = readh5_imputation(paste0(ds_dir, "/VIPER.hdf5"))

[1] "./data/MELANOMA/ds_0.5/r1/VIPER.hdf5"
[1] 12509  8498

data_list[["scImpute"]] = readh5_imputation(paste0(ds_dir, "/scImpute.hdf5"))

[1] "./data/MELANOMA/ds_0.5/r1/scImpute.hdf5"
[1] 12509  8498

cell_number = ncol(data_list[["Raw"]])
gene_number = length(used_genes)

Settings

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

MAE

MAE_mat = matrix(nrow = cell_number, ncol = length(method_names), dimnames = list(c(), method_names))
ls_raw = colSums(data_list[["Raw"]])
for(ii in method_names){
  ls_this = colSums(data_list[[ii]])
    scale_factor = ls_raw / ls_this
    MAE_cell = sapply(seq(cell_number), function(x){
      expressed_mask = data_list[["Raw"]][, x] > 0
      expressed_number = sum(expressed_mask)
      error = data_list[["Raw"]][expressed_mask, x] - (data_list[[ii]][expressed_mask, x] * scale_factor[x])
      return(sum(abs(error)) / expressed_number)
    })
  MAE_mat[, ii] = MAE_cell
  print(ii)
}

[1] "Observed"
[1] "DISC"
[1] "scVI"
[1] "MAGIC"
[1] "DCA"
[1] "scScope"
[1] "DeepImpute"
[1] "VIPER"
[1] "scImpute"

MAE_mat = MAE_mat[rowSums(is.na(MAE_mat)) < 1, ]
MAE_df = melt(t(MAE_mat))
MAE_levels = colnames(MAE_mat)[c(1, order(colMeans(MAE_mat)[-1], decreasing = F) + 1)]

ggplot(MAE_df, aes(x = factor(Var1, levels = MAE_levels), y = value, fill = factor(Var1, levels = MAE_levels))) +
  geom_boxplot(outlier.shape = NA) + stat_boxplot(geom = "errorbar", width = 0.3) +
  ylim(min(apply(MAE_mat, 2, quantile, 0.1)), max(apply(MAE_mat, 2, quantile, 0.9))) + theme_classic() +
  labs(x="", y="MAE") + scale_fill_manual(values=method_color) +
  theme(axis.text.x = element_text(size = 12,angle = 45, hjust = 1, vjust = 1, face = "bold"),
        axis.text.y = element_text(size = 12, hjust = 1, vjust = 1, face = "bold"),
        axis.title.y = element_text(size = 15, face = "bold"),
        legend.position = "none")

CMD

vst_file = paste0(output_dir, "/vst_gene.tsv")
if(file.exists(vst_file)){
  hvg_info = read.table(vst_file)
  print("load vst_file")
}else{
  hvg_info = FindVariableFeatures_vst_by_genes(data_list[["Raw"]])
  hvg_info = hvg_info[order(hvg_info$variance.standardized, decreasing = T), ]
  write.table(hvg_info, paste0(output_dir, "/vst_gene.tsv"), sep = "\t", quote = F, row.names = T, col.names = T)
}
used_feature_genes = rownames(hvg_info)[1:300]
cor_all = list()
for(ii in names(data_list)){
  if(ii == "Raw"){
    cor_all[[ii]] = calc_cor_mat(data_list[[ii]][used_feature_genes, ])
  }else{
    cor_all[[ii]] = calc_cor_mat(delete_lt0.5(data_list[[ii]])[used_feature_genes, ])
  }
  print(ii)
}

[1] "Raw"
[1] "Observed"
[1] "DISC"
[1] "scVI"
[1] "MAGIC"
[1] "DCA"
[1] "scScope"
[1] "DeepImpute"
[1] "VIPER"
[1] "scImpute"

saveRDS(cor_all, paste0(output_dir, "/cor_all.rds"))
cmd_value = c()
for(ii in method_names){
  cmd_value = c(cmd_value, calc_cmd(cor_all[["Raw"]], cor_all[[ii]]))
}
names(cmd_value) = method_names

barplot_usage(cmd_value, main = "", cex.main = 1.5, bar_color = method_color, text_color = text_color, use_data_order = T, ylab = "CMD", cex.lab = 1.5, font.main = 1, ylim = c(-0.1, 1), use_border = F)

Gene correlation

gene_corr_mat = matrix(nrow = gene_number, ncol = length(method_names), dimnames = list(c(), method_names))
for(ii in method_names){
  gene_corr_mat[, ii] = calc_corr(data_list[["Raw"]], data_list[[ii]], "gene")
  print(ii)
}

[1] "Observed"
[1] "DISC"
[1] "scVI"
[1] "MAGIC"
[1] "DCA"
[1] "scScope"
[1] "DeepImpute"
[1] "VIPER"
[1] "scImpute"

gene_corr_mat = gene_corr_mat[rowSums(is.na(gene_corr_mat)) < 1, ]
gene_corr_df = melt(t(gene_corr_mat))
gene_corr_levels = colnames(gene_corr_mat)[c(1, order(colMeans(gene_corr_mat)[-1], decreasing = T) + 1)]

ggplot(gene_corr_df, aes(x = factor(Var1, levels = gene_corr_levels), y = value, fill = factor(Var1, levels = gene_corr_levels))) +
  geom_boxplot(outlier.shape = NA) + stat_boxplot(geom = "errorbar", width = 0.3) +
  ylim(min(apply(gene_corr_mat, 2, quantile, 0.1)), max(apply(gene_corr_mat, 2, quantile, 0.9))) + theme_classic() +
  labs(x="", y="Gene correlation with reference") + scale_fill_manual(values=method_color) +
  theme(axis.text.x = element_text(size = 12,angle = 45, hjust = 1, vjust = 1, face = "bold"),
        axis.text.y = element_text(size = 12, hjust = 1, vjust = 1, face = "bold"),
        axis.title.y = element_text(size = 15, face = "bold"),
        legend.position = "none")

Cell correlation

cell_corr_mat = matrix(nrow = cell_number, ncol = length(method_names), dimnames = list(c(), method_names))
for(ii in method_names){
  cell_corr_mat[, ii] = calc_corr(data_list[["Raw"]], data_list[[ii]], "cell")
  print(ii)
}

[1] "Observed"
[1] "DISC"
[1] "scVI"
[1] "MAGIC"
[1] "DCA"
[1] "scScope"
[1] "DeepImpute"
[1] "VIPER"
[1] "scImpute"

cell_corr_mat = cell_corr_mat[rowSums(is.na(cell_corr_mat)) < 1, ]
cell_corr_df = melt(t(cell_corr_mat))
cell_corr_levels = colnames(cell_corr_mat)[c(1, order(colMeans(cell_corr_mat)[-1], decreasing = T) + 1)]

ggplot(cell_corr_df, aes(x = factor(Var1, levels = cell_corr_levels), y = value, fill = factor(Var1, levels = cell_corr_levels))) +
  geom_boxplot(outlier.shape = NA) + stat_boxplot(geom = "errorbar", width = 0.3) +
  ylim(min(apply(cell_corr_mat, 2, quantile, 0.1)), max(apply(cell_corr_mat, 2, quantile, 0.9))) + theme_classic() +
  labs(x="", y="Cell correlation with reference") + scale_fill_manual(values=method_color) +
  theme(axis.text.x = element_text(size = 12,angle = 45, hjust = 1, vjust = 1, face = "bold"),
        axis.text.y = element_text(size = 12, hjust = 1, vjust = 1, face = "bold"),
        axis.title.y = element_text(size = 15, face = "bold"),
        legend.position = "none")