Hierarchical annotation of immune cells in scRNA-Seq data based on ssGSEA algorithm. Fork for large datasets with QOL improvements.
byarielm.fyi
1#' @title Seurat Object.
2#' @description Creating dynamic Seurat Object to add annotation information
3#' and draw images.
4#' @details Input takes a count, genelist and scImmuCC result, returns tsne,
5#' umap, Dotplot and pheatmaps.
6#' @param count a matrix with cell unique barcodes as column names and
7#' gene names as row names
8#' @param genematrix Data frame with cell types as column names
9#' @param ssGSEA_result Data frame, scImmuCC return result
10#' @param filename Custom file name, char
11#' @return 4 pictures
12#' @import ggplot2
13#' @import Seurat
14#' @importFrom Seurat CreateSeuratObject NormalizeData FindVariableFeatures
15#' ScaleData PercentageFeatureSet VariableFeatures RunPCA RunUMAP RunTSNE
16#' DimPlot DotPlot DoHeatmap ProjectDim JackStraw ScoreJackStraw
17
18seurat_Heatmap <- function(count, genematrix, ssGSEA_result, filename) {
19 count <- count[, !duplicated(colnames(count))]
20 celltype <- intersect(colnames(count), ssGSEA_result[, 1])
21 labels <- ssGSEA_result[which(ssGSEA_result[, 1] %in% celltype), ]
22 count <- count[, labels[, 1]]
23
24 seurat.data <- CreateSeuratObject(counts = count, project = filename)
25 seurat.data
26 seurat.data[["percent.mt"]] <- PercentageFeatureSet(seurat.data,
27 pattern = "^MT-"
28 )
29 HB.genes_total <- c(
30 "HBA1", "HBA2", "HBB", "HBD", "HBE1", "HBG1", "HBG2",
31 "HBM", "HBQ1", "HBZ"
32 )
33 HB_m <- match(HB.genes_total, rownames(seurat.data@assays$RNA))
34
35 HB.genes <- rownames(seurat.data@assays$RNA)[HB_m]
36 HB.genes <- HB.genes[!is.na(HB.genes)]
37 seurat.data[["percent.HB"]] <- PercentageFeatureSet(seurat.data,
38 features = HB.genes
39 )
40
41 seurat.data <- NormalizeData(seurat.data,
42 normalization.method = "LogNormalize",
43 scale.factor = 10000
44 )
45 seurat.data <- FindVariableFeatures(seurat.data,
46 selection.method = "vst",
47 nfeatures = 2000
48 )
49 all.genes <- rownames(seurat.data)
50 seurat.data <- ScaleData(seurat.data, features = all.genes)
51
52 # Perform linear dimensional reduction
53 counts <- seurat.data[["RNA"]]$counts
54 cells <- length(counts[2, ])
55 if (cells > 50) {
56 seurat.data <- RunPCA(seurat.data,
57 features = VariableFeatures(object = seurat.data)
58 )
59 } else {
60 seurat.data <- RunPCA(seurat.data,
61 npcs = (cells - 1),
62 features = VariableFeatures(object = seurat.data)
63 )
64 }
65
66 head(seurat.data@reductions$pca@cell.embeddings)
67 head(seurat.data@reductions$pca@feature.loadings)
68 seurat.data <- ProjectDim(object = seurat.data)
69 seurat.data <- JackStraw(seurat.data, num.replicate = 100)
70 seurat.data <- ScoreJackStraw(seurat.data, dims = 1:20)
71
72 seurat.data <- RunUMAP(seurat.data, dims = 1:10)
73 n <- length(count[2, ])
74
75 min_cell_count <- 100
76 if (ncol(seurat.data) >= min_cell_count) {
77 seurat.data <- RunTSNE(seurat.data, dims = 1:10, check_duplicates = FALSE)
78 } else {
79 seurat.data <- RunTSNE(seurat.data,
80 dims = 1:10, perplexity = 5,
81 check_duplicates = FALSE
82 )
83 }
84
85 head(seurat.data@reductions$tsne@cell.embeddings)
86
87 # Add annotation information to Seurat object
88 seurat.data@meta.data$cell_type_pred <- labels[, 2]
89
90 pdf(paste(filename, "_tSNE", ".pdf", sep = ""), width = 12, height = 10)
91 p1 <- DimPlot(seurat.data,
92 reduction = "tsne",
93 group.by = "cell_type_pred", label = TRUE, pt.size = 1
94 )
95 plot(p1)
96 dev.off()
97 pdf(paste(filename, "_UMAP", ".pdf", sep = ""), width = 12, height = 10)
98 p2 <- DimPlot(seurat.data,
99 reduction = "umap",
100 group.by = "cell_type_pred", label = TRUE, pt.size = 1
101 )
102 plot(p2)
103 dev.off()
104
105 genelist <- as.list(genematrix)
106 genelist <- lapply(genelist, function(x) x[!is.na(x)])
107 gene <- c()
108 for (i in 1:length(genelist)) {
109 gene <- c(gene, genelist[[i]])
110 }
111
112 length(gene)
113 gene <- unique(gene)
114 Y <- intersect(gene, rownames(count))
115 length(Y)
116 features <- Y
117 pdf(paste(filename, "_UMAP_DotPlot", ".pdf", sep = ""),
118 width = 12,
119 height = 10
120 )
121 p3 <- DotPlot(seurat.data,
122 features = features,
123 group.by = "cell_type_pred"
124 ) + RotatedAxis()
125 plot(p3)
126 dev.off()
127
128 pdf(paste(filename, "_DoHeatmap.pdf", sep = ""), width = 12, height = 10)
129 p4 <- DoHeatmap(seurat.data,
130 features = Y,
131 group.by = "cell_type_pred"
132 ) + NoLegend()
133 plot(p4)
134 dev.off()
135}