Data preparation for 10X_5CL

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.

Generate data file

The original 10X_5CL data can be found at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE126906.
The 10X_5CL metadata can be found at https://github.com/LuyiTian/sc_mixology/blob/master/data/csv/sc_10x_5cl.metadata.csv.gz.
We use 3918 filtered cells and all genes to make raw.loom (genes will be filtered in the following step gene selection as our DISC paper).

if(!file.exists("./data/10X_5CL/raw.loom")){
  original_data = as.matrix(read.csv("./data/10X_5CL/original_data/GSM3618014_gene_count.csv.gz", header = T, sep = ",", row.names = 1))# 32895  5001
  gz_path = "./data/annotation/Homo_sapiens.GRCh37.87.gtf.gz"
  annotation_mat = get_map(gz_path)
  mapping = annotation_mat$gene_name
  names(mapping) = annotation_mat$gene_id
  gene_name_candidate = mapping[rownames(original_data)]
  gene_name = rownames(original_data)
  gene_name[!is.na(gene_name_candidate)] = gene_name_candidate[!is.na(gene_name_candidate)]
  for(ii in gene_name[duplicated(gene_name)]){
    index = 0
    for(jj in which(gene_name == ii)){
      gene_name[jj] = paste0(gene_name[jj], "_duplicate_", index)
      index = index + 1
    }
  }
  rownames(original_data) = gene_name
  print(dim(original_data))
  save_h5("./data/10X_5CL/original.loom", t(original_data))
  metadata = as.matrix(read.csv("./data/10X_5CL/original_data/sc_10x_5cl.metadata.csv.gz", header = T, row.names = 1))
  cell_id = rownames(metadata)
  raw_data = as.matrix(original_data[, cell_id])
  cell_id_new = paste(metadata[, "cell_line"], cell_id, sep = "_")
  colnames(raw_data) = cell_id_new
  print(dim(raw_data))
  save_h5("./data/10X_5CL/raw.loom", t(raw_data))
}else{
  raw_data = readh5_loom("./data/10X_5CL/raw.loom")
  cell_id = colnames(raw_data)
  gene_name = rownames(raw_data)
  print(dim(raw_data))
}

Treat input as file
|--------------------------------------------------|
|==================================================|

[1] 32895  5001
[1] 32895  3918

print("10X_5CL...OK!")

[1] "10X_5CL...OK!"

The bulk data can be found at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE86337.

if(!file.exists("./data/10X_5CL/bulk.loom")){
  library(org.Hs.eg.db)
  gene_bulk_mat_raw = as.matrix(read.table("./data/10X_5CL/original_data/GSE86337_reverse.stranded.unfiltered.count.matrix.txt.gz", header = T, row.names = 1, sep = "\t"))
  sample_mapping = as.matrix(read.csv("./data/10X_5CL/original_data/GSE86337_sample_id_sample_name.txt", sep = "\t"))
  sample_name = sample_mapping[, 2]
  names(sample_name) = sample_mapping[, 1]
  colnames(gene_bulk_mat_raw) = sample_name[colnames(gene_bulk_mat_raw)]
  gene_name_all = select(org.Hs.eg.db, key = as.character(rownames(gene_bulk_mat_raw)), columns = c("SYMBOL"), keytype = "ENTREZID")$SYMBOL
  rownames(gene_bulk_mat_raw) = gene_name_all
  gene_bulk_filt_0 = gene_bulk_mat_raw[!is.na(gene_name_all), ]
  gene_name_used = rownames(gene_bulk_filt_0)
  used_index = sapply(unique(gene_name_used), function(x){
    this_index = which(gene_name_used == x)
    if(length(this_index) > 1){
      return(this_index[which.max(rowSums(gene_bulk_filt_0[this_index, ]))])
    }else{
      return(this_index)
    }
  })
  gene_bulk_mat = gene_bulk_filt_0[used_index, ]
  save_h5("./data/10X_5CL/bulk.loom", as.matrix(t(gene_bulk_mat)))
}else{
  gene_bulk_mat = readh5_loom("./data/10X_5CL/bulk.loom")
}

[1] "./data/10X_5CL/bulk.loom"
[1] 23753    10

print(dim(gene_bulk_mat))

[1] 23753    10

print("10X_5CL_Bulk...OK!")

[1] "10X_5CL_Bulk...OK!"