createSCE_kallisto.R

#
# Load libraries
#
library(magrittr)
library(data.table)
library(Matrix)
library(biomaRt)
library(SingleCellExperiment)

#
# Read tx2gene file
#
tx2gene <- read.table(snakemake@input$tx2gene, sep="\t", col.names = c("transcript", "gene"))

#
# Read quant file from Kallisto and construct expression matrix
#
sparse_data <- read.table(snakemake@input$abundance, header = F, sep = "\t")
umi_counts <- Matrix::sparseMatrix(i = sparse_data$V1+1, j = sparse_data$V2+1, x = sparse_data$V3)

cells <- readLines(snakemake@input$cells)
colnames(umi_counts) <- cells
rownames(umi_counts) <- paste0("EC", 1:nrow(umi_counts))

#
# Provide mapping between ECs and transcripts/genes
#
data.table::fread(snakemake@input$ec, data.table = F) %>%
  magrittr::extract2(2) %>%
  strsplit(",", fixed = T) %>%
  lapply(as.numeric) %>%
  lapply(magrittr::add, 1) -> ec_members

#
# Assemble list of expression matrices and adjust colnames
#
expression <- list(umi = umi_counts)

#
# Annotate gene names
#
if (snakemake@config$container$annotate) {
  ec_members %>%
    lapply(function(m) {unique(tx2gene[m, "gene"])}) -> ec2gene

  bm <- biomaRt::getBM(c(snakemake@config$container$biomart$filter, 'external_gene_name'),
                       filters = snakemake@config$container$biomart$filter,
                       values = ec2gene %>% unlist %>% unique,
                       mart = biomaRt::useMart("ensembl", dataset = snakemake@config$container$biomart$dataset))

  gene2symbol <- bm[, 'external_gene_name']
  names(gene2symbol) <- bm[, snakemake@config$container$biomart$filter]

  rowData_ <- data.frame(transcript_ids = sapply(ec_members, function(m) paste0(tx2gene[m, "transcript"], collapse = ",")),
                         gene_ids = sapply(ec2gene, function(m) paste0(m, collapse = ",")),
                         symbols = sapply(ec2gene, function(m) { paste0(unname(gene2symbol[m]), collapse = ",") }))
} else {
  rowData_ <- NULL
}

#
# Assemble metadata
#
metadata_files <- snakemake@input$metadata
if (!is.null(metadata_files)) {
  metadata_ <- lapply(metadata_files, function(p) metadata <- data.table::fread(input = p, sep = "\t", header = T, stringsAsFactors = T, na.strings = "NA"))
  metadata <- Reduce(function(a, b) dplyr::left_join(a, b, by = snakemake@config$samplesheet$index), metadata_)

  m <- match(colnames(expression[[1]]), make.names(metadata[[which(colnames(metadata) == snakemake@config$samplesheet$index)]]))
  colData_ <- as.data.frame(metadata[na.omit(m), ])
  rownames(colData_) <- make.names(colData_[, snakemake@config$samplesheet$index])
} else {
  colData_ <- NULL
}

if (!is.null(colData_) & !is.null(rowData_)) {
  scData <- SingleCellExperiment::SingleCellExperiment(assays = expression, colData = colData_, rowData = rowData_)
}
if (is.null(colData_) & is.null(rowData_)) {
  scData <- SingleCellExperiment::SingleCellExperiment(assays = expression)
}
if (is.null(colData_) & !is.null(rowData_)) {
  scData <- SingleCellExperiment::SingleCellExperiment(assays = expression, rowData = rowData_)
}
if (is.null(rowData_) & !is.null(colData_)) {
  scData <- SingleCellExperiment::SingleCellExperiment(assays = expression, colData = colData_)
}

#
# Save generated object
#
save(scData, file = snakemake@output$rdata)