evaluation.R

require(Seurat)
require(scmap)
require(DropletUtils)
require(SeuratConverter)

initLung <- function() {
  cache <- '../dataset 1/SC01+02.Robj'
  if (!file.exists(cache)) {
    sc01.data <- Read10X(data.dir = "../dataset 1/SC01/")
    sc01 <- CreateSeuratObject(raw.data = sc01.data, min.cells = 3, min.genes = 200, project = "SC01")
    [email protected] <- paste("SC01", [email protected], sep = "_")
    colnames(x = [email protected]) <- paste("SC01", colnames(x = [email protected]), sep = "_")
    rownames(x = [email protected]) <- paste("SC01", rownames(x = [email protected]), sep = "_")
    
    sc02.data <- Read10X(data.dir = "../dataset 1/SC02/")
    lung <- AddSamples(object = sc01, new.data = sc02.data, add.cell.id = "SC02")
    
    mito.genes <- grep(pattern = "^mt-", x = rownames(x = lung@data), value = TRUE)
    percent.mito <- Matrix::colSums([email protected][mito.genes, ])/Matrix::colSums([email protected])
    lung <- AddMetaData(object = lung, metadata = percent.mito, col.name = "percent.mito")
    lung <- FilterCells(object = lung, subset.names = c("nGene", "percent.mito"), 
                        low.thresholds = c(300, -Inf), high.thresholds = c(4000, 0.1))
    lung <- NormalizeData(object = lung, normalization.method = "LogNormalize", scale.factor = 10000)
    lung <- ScaleData(object = lung, vars.to.regress = c("nUMI", "percent.mito"))
    lung <- FindVariableGenes(object = lung, mean.function = ExpMean, dispersion.function = LogVMR, 
                              x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)
    lung <- RunPCA(object = lung, pc.genes = [email protected], do.print = TRUE, pcs.print = 1:5, 
                   genes.print = 5, pcs.compute = 40)
    lung <- ProjectPCA(object = lung, do.print = FALSE)
    lung <- FindClusters(object = lung, reduction.type = "pca", dims.use = 1:35, 
                         resolution = 0.5, print.output = 0, save.SNN = TRUE, force.recalc = T)
    lung <- RunTSNE(object = lung, dims.use = 1:35, do.fast = TRUE, check_duplicates = FALSE)
    
    save(lung, file=cache)
  } else {
    load(cache)
  }
  return(lung)
}

initSC03 <- function() {
  cache <- '../dataset 1/sc03.Robj'
  if (!file.exists(cache)) {
    sc03.data <- Read10X(data.dir = "../dataset 1/SC03/")
    sc03.seurat <- CreateSeuratObject(raw.data = sc03.data, min.cells = 3, min.genes = 200, project = "SC03")
    sc03.mito.genes <- grep(pattern = "^mt-", x = rownames(x = sc03.seurat@data), value = TRUE)
    sc03.percent.mito <- Matrix::colSums([email protected][sc03.mito.genes, ])/Matrix::colSums([email protected])
    sc03.seurat <- AddMetaData(object = sc03.seurat, metadata = sc03.percent.mito, col.name = "percent.mito")
    sc03.seurat <- FilterCells(object = sc03.seurat, subset.names = c("nGene", "percent.mito"), 
                        low.thresholds = c(300, -Inf), high.thresholds = c(4000, 0.1))
    sc03.seurat <- NormalizeData(object = sc03.seurat, normalization.method = "LogNormalize", scale.factor = 10000)
    sc03.seurat <- ScaleData(object = sc03.seurat, vars.to.regress = c("nUMI", "percent.mito"))
    sc03.seurat <- FindVariableGenes(object = sc03.seurat, mean.function = ExpMean, dispersion.function = LogVMR, 
                              x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)
    sc03.seurat <- RunPCA(object = sc03.seurat, pc.genes = [email protected], do.print = TRUE, pcs.print = 1:5, 
                   genes.print = 5, pcs.compute = 40)
    sc03.seurat <- ProjectPCA(object = sc03.seurat, do.print = FALSE)
    sc03.seurat <- FindClusters(object = sc03.seurat, reduction.type = "pca", dims.use = 1:21, 
                         resolution = 0.5, print.output = 0, save.SNN = TRUE, force.recalc = T)
    sc03.seurat <- RunTSNE(object = sc03.seurat, dims.use = 1:21, do.fast = TRUE, check_duplicates = FALSE)
    save(sc03.seurat, file=cache)
  } else {
    load(cache)
  }
  return(sc03.seurat)
}

initSC01 <- function(lung) {
  cache <- '../dataset 1/sc01-sce.Robj'
  if (!file.exists(cache)) {
    sc01 <- as(SubsetData(lung, [email protected][[email protected]$orig.ident=="SC01"]), "SingleCellExperiment")
    counts(sc01) <- as.matrix(assay(sc01, "raw.data"))
    logcounts(sc01) <- log2(counts(sc01) + 1)
    rowData(sc01)$feature_symbol <- rownames(sc01)
    colData(sc01)$cell_type1 <- lung@ident[[email protected]$orig.ident == "SC01"]
    sc01 <- selectFeatures(sc01)
    sc01 <- indexCluster(sc01)
    sc01 <- indexCell(sc01)
    saveRDS(sc01, file=cache)
    return(sc01)
  } else {
    return(readRDS(cache))
  }
}

initSC03sce <- function(sc03.seurat) {
  cache <- '../dataset 1/sc03-sce.Robj'
  if (!file.exists(cache)) {
    sc03 <- as(sc03.seurat, "SingleCellExperiment")
    counts(sc03) <- as.matrix(assay(sc03, "raw.data"))
    logcounts(sc03) <- log2(counts(sc03) + 1)
    rowData(sc03)$feature_symbol <- rownames(sc03)
    colData(sc03)$cell_type1 <- sc03.seurat@ident
    saveRDS(sc03, file=cache)
    return(sc03)
  } else {
    return(readRDS(cache))
  }
}

runCluster <- function(sc03, sc01) {
  scmapCluster_results <- scmapCluster(projection = sc03, 
                                       index_list = list(sc01=metadata(sc01)$scmap_cluster_index))
  plot(getSankey(colData(sc03)$cell_type1, scmapCluster_results$combined_labs))
}

runCell <- function(sc03, sc01) {
  scmapCell_results <- scmapCell(projection=sc03, 
                                 index_list=list(sc01=metadata(sc01)$scmap_cell_index))
  scmapCell_clusters <- scmapCell2Cluster(scmapCell_results, list(colData(sc01)$cell_type1))
  plot(getSankey(colData(sc03)$cell_type1, scmapCell_clusters$combined_labs))
}

lung <- initLung()
sc01 <- initSC01(lung)
sc03s <- initSC03()
sc03 <- initSC03sce(sc03s)
runCluster(sc03, sc01)
runCell(sc03, sc01)