createSCE_salmon.R

#
# Load libraries
#
library(tximport)
library(biomaRt)
library(SingleCellExperiment)

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

#
# Read quant files from Salmon
#
quant_files <- snakemake@input$quant_files
names(quant_files) <- make.names(snakemake@params$samples)

txi_salmon <- tximport::tximport(quant_files, type = "salmon", tx2gene = tx2gene)

#
# Assemble list of expression matrices and adjust colnames
#
expression <- txi_salmon[c("abundance", "counts")]

#
# Exclude genes that show zero expression across all cells
#
if (snakemake@config$container$exclude_zeros) {
  keep <- rowSums(expression$counts) > 0
  expression <- lapply(expression, function(e) e[keep,])
}

#
# Annotate gene names
#
if (snakemake@config$container$annotate) {
  identifier <- rownames(expression[[1]])
  bm <- biomaRt::getBM(snakemake@config$container$biomart$fields,
                       filters = snakemake@config$container$biomart$filter,
                       values = identifier,
                       mart = biomaRt::useMart("ensembl", dataset = snakemake@config$container$biomart$dataset))
  m <- match(identifier, bm[, snakemake@config$container$biomart$filter])

  data <- data.frame(bm[na.omit(m), ])
  rownames(data) <- data[, snakemake@config$container$biomart$filter]
  colnames(data) <- snakemake@config$container$biomart$fields

  rowData_ <- data[,c(2:ncol(data))]

  if (nrow(rowData_) != nrow(expression[[1]])) {
    warning("Could not annotate all features using biomaRt. Omitting unannotated features.")
    expression <- lapply(expression, function(l) l[which(!is.na(m)), ])
  }
} 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[, 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)