Mercurial > repos > iuc > seurat

comparison Seurat.R @ 6:764f076e9d52 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/seurat commit e8eeff2efea68e1e5bc697d0ff6a5808cd978db2"
author iuc
date Thu, 23 Jul 2020 11:22:05 -0400
parents 321bdd834266
children c4db6ec33fec
comparison
equal deleted inserted replaced
5:06ed31cf52ed 6:764f076e9d52
8 #' low_thresholds: "" 8 #' low_thresholds: ""
9 #' high_thresholds: "" 9 #' high_thresholds: ""
10 #' numPCs: "" 10 #' numPCs: ""
11 #' cells_use: "" 11 #' cells_use: ""
12 #' resolution: "" 12 #' resolution: ""
13 #' perplexity: ""
13 #' min_pct: "" 14 #' min_pct: ""
14 #' logfc_threshold: "" 15 #' logfc_threshold: ""
15 #' showcode: "" 16 #' showcode: ""
16 #' warn: "" 17 #' warn: ""
17 #' varstate: "" 18 #' varstate: ""
20 #' PCplots: "" 21 #' PCplots: ""
21 #' tsne: "" 22 #' tsne: ""
22 #' heatmaps: "" 23 #' heatmaps: ""
23 #' --- 24 #' ---
24 25
26 # nolint start
25 #+ echo=F, warning = F, message=F 27 #+ echo=F, warning = F, message=F
26 options(show.error.messages = F, error = function(){cat(geterrmessage(), file = stderr()); q("no", 1, F)}) 28 options(show.error.messages = F, error = function() {
29 cat(geterrmessage(), file = stderr()); q("no", 1, F)
30 })
27 showcode <- as.logical(params$showcode) 31 showcode <- as.logical(params$showcode)
28 warn <- as.logical(params$warn) 32 warn <- as.logical(params$warn)
29 varstate <- as.logical(params$varstate) 33 varstate <- as.logical(params$varstate)
30 vlnfeat <- as.logical(params$vlnfeat) 34 vlnfeat <- as.logical(params$vlnfeat)
31 featplot <- as.logical(params$featplot) 35 featplot <- as.logical(params$featplot)
32 PCplots <- as.logical(params$PCplots) 36 pc_plots <- as.logical(params$PCplots)
33 tsne <- as.logical(params$tsne) 37 tsne <- as.logical(params$tsne)
34 heatmaps <- as.logical(params$heatmaps) 38 heatmaps <- as.logical(params$heatmaps)
35 39
36 # we need that to not crash Galaxy with an UTF-8 error on German LC settings. 40 # we need that to not crash Galaxy with an UTF-8 error on German LC settings.
37 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") 41 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
38 42
39
40 #+ echo = F, warning = `warn`, include =`varstate` 43 #+ echo = F, warning = `warn`, include =`varstate`
41 min_cells <- as.integer(params$min_cells) 44 min_cells <- as.integer(params$min_cells)
42 min_genes <- as.integer(params$min_genes) 45 min_genes <- as.integer(params$min_genes)
43 low_thresholds <- as.integer(params$low_thresholds) 46 low_thresholds <- as.integer(params$low_thresholds)
44 high_thresholds <- as.integer(params$high_thresholds) 47 high_thresholds <- as.integer(params$high_thresholds)
45 numPCs <- as.integer(params$numPCs) 48 num_pcs <- as.integer(params$numPCs)
46 cells_use <- as.integer(params$cells_use) 49 cells_use <- as.integer(params$cells_use)
47 resolution <- as.double(params$resolution) 50 resolution <- as.double(params$resolution)
51 perplexity <- as.integer(params$perplexity)
48 min_pct <- as.double(params$min_pct) 52 min_pct <- as.double(params$min_pct)
49 logfc_threshold <- as.double(params$logfc_thresh) 53 logfc_threshold <- as.double(params$logfc_thresh)
50 print(paste0("Minimum cells: ", min_cells)) 54 print(paste0("Minimum cells: ", min_cells))
51 print(paste0("Minimum features: ", min_genes)) 55 print(paste0("Minimum features: ", min_genes))
52 print(paste0("Umi low threshold: ", low_thresholds)) 56 print(paste0("Umi low threshold: ", low_thresholds))
53 print(paste0("Umi high threshold: ", high_thresholds)) 57 print(paste0("Umi high threshold: ", high_thresholds))
54 print(paste0("Number of principal components: ", numPCs)) 58 print(paste0("Number of principal components: ", num_pcs))
55 print(paste0("Resolution: ", resolution)) 59 print(paste0("Resolution: ", resolution))
60 print(paste0("Perplexity: ", perplexity))
56 print(paste0("Minimum percent of cells", min_pct)) 61 print(paste0("Minimum percent of cells", min_pct))
57 print(paste0("Logfold change threshold", logfc_threshold)) 62 print(paste0("Logfold change threshold", logfc_threshold))
58 63
59 #+ echo = FALSE 64 #+ echo = FALSE
60 if(showcode == TRUE) print("Read in data, generate inital Seurat object") 65 if (showcode == TRUE) print("Read in data, generate inital Seurat object")
61 #+ echo = `showcode`, warning = `warn`, message = F 66 #+ echo = `showcode`, warning = `warn`, message = F
62 counts <- read.delim(params$counts, row.names = 1) 67 counts <- read.delim(params$counts, row.names = 1)
63 seuset <- Seurat::CreateSeuratObject(counts = counts, min.cells = min_cells, min.features = min_genes) 68 seuset <- Seurat::CreateSeuratObject(counts = counts, min.cells = min_cells, min.features = min_genes)
64 69
65 #+ echo = FALSE 70 #+ echo = FALSE
66 if(showcode == TRUE && vlnfeat == TRUE) print("Raw data vizualization") 71 if (showcode == TRUE && vlnfeat == TRUE) print("Raw data vizualization")
67 #+ echo = `showcode`, warning = `warn`, include=`vlnfeat` 72 #+ echo = `showcode`, warning = `warn`, include=`vlnfeat`
68 Seurat::VlnPlot(object = seuset, features = c("nFeature_RNA", "nCount_RNA")) 73 Seurat::VlnPlot(object = seuset, features = c("nFeature_RNA", "nCount_RNA"))
69 Seurat::FeatureScatter(object = seuset, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") 74 Seurat::FeatureScatter(object = seuset, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
70 75
71 #+ echo = FALSE 76 #+ echo = FALSE
72 if(showcode == TRUE) print("Filter and normalize for UMI counts") 77 if (showcode == TRUE) print("Filter and normalize for UMI counts")
73 #+ echo = `showcode`, warning = `warn` 78 #+ echo = `showcode`, warning = `warn`
74 seuset <- subset(seuset, subset = `nCount_RNA` > low_thresholds & `nCount_RNA` < high_thresholds) 79 seuset <- subset(seuset, subset = `nCount_RNA` > low_thresholds & `nCount_RNA` < high_thresholds)
75 seuset <- Seurat::NormalizeData(seuset, normalization.method = "LogNormalize", scale.factor = 10000) 80 seuset <- Seurat::NormalizeData(seuset, normalization.method = "LogNormalize", scale.factor = 10000)
76 81
77 #+ echo = FALSE 82 #+ echo = FALSE
78 if(showcode == TRUE && featplot == TRUE) print("Variable Genes") 83 if (showcode == TRUE && featplot == TRUE) print("Variable Genes")
79 #+ echo = `showcode`, warning = `warn`, include = `featplot` 84 #+ echo = `showcode`, warning = `warn`, include = `featplot`
80 seuset <- Seurat::FindVariableFeatures(object = seuset, selection.method = "mvp") 85 seuset <- Seurat::FindVariableFeatures(object = seuset, selection.method = "mvp")
81 Seurat::VariableFeaturePlot(seuset, cols = c("black", "red"), selection.method = "disp") 86 Seurat::VariableFeaturePlot(seuset, cols = c("black", "red"), selection.method = "disp")
82 seuset <- Seurat::ScaleData(object = seuset, vars.to.regress = "nCount_RNA") 87 seuset <- Seurat::ScaleData(object = seuset, vars.to.regress = "nCount_RNA")
83 88
84 #+ echo = FALSE 89 #+ echo = FALSE
85 if(showcode == TRUE && PCplots == TRUE) print("PCA Visualization") 90 if (showcode == TRUE && pc_plots == TRUE) print("PCA Visualization")
86 #+ echo = `showcode`, warning = `warn`, include = `PCplots` 91 #+ echo = `showcode`, warning = `warn`, include = `pc_plots`
87 seuset <- Seurat::RunPCA(seuset, npcs = numPCs) 92 seuset <- Seurat::RunPCA(seuset, npcs = num_pcs)
88 Seurat::VizDimLoadings(seuset, dims = 1:2) 93 Seurat::VizDimLoadings(seuset, dims = 1:2)
89 Seurat::DimPlot(seuset, dims = c(1,2), reduction = "pca") 94 Seurat::DimPlot(seuset, dims = c(1, 2), reduction = "pca")
90 Seurat::DimHeatmap(seuset, dims = 1:numPCs, nfeatures = 30, reduction = "pca") 95 Seurat::DimHeatmap(seuset, dims = 1:num_pcs, nfeatures = 30, reduction = "pca")
91 seuset <- Seurat::JackStraw(seuset, dims=numPCs, reduction = "pca", num.replicate = 100) 96 seuset <- Seurat::JackStraw(seuset, dims = num_pcs, reduction = "pca", num.replicate = 100)
92 seuset <- Seurat::ScoreJackStraw(seuset, dims = 1:numPCs) 97 seuset <- Seurat::ScoreJackStraw(seuset, dims = 1:num_pcs)
93 Seurat::JackStrawPlot(seuset, dims = 1:numPCs) 98 Seurat::JackStrawPlot(seuset, dims = 1:num_pcs)
94 Seurat::ElbowPlot(seuset, ndims = numPCs, reduction = "pca") 99 Seurat::ElbowPlot(seuset, ndims = num_pcs, reduction = "pca")
95 100
96 #+ echo = FALSE 101 #+ echo = FALSE
97 if(showcode == TRUE && tsne == TRUE) print("tSNE") 102 if (showcode == TRUE && tsne == TRUE) print("tSNE")
98 #+ echo = `showcode`, warning = `warn`, include = `tsne` 103 #+ echo = `showcode`, warning = `warn`, include = `tsne`
99 seuset <- Seurat::FindNeighbors(object = seuset) 104 seuset <- Seurat::FindNeighbors(object = seuset)
100 seuset <- Seurat::FindClusters(object = seuset) 105 seuset <- Seurat::FindClusters(object = seuset)
101 seuset <- Seurat::RunTSNE(seuset, dims = 1:numPCs, resolution = resolution) 106 if (perplexity == -1) {
107 seuset <- Seurat::RunTSNE(seuset, dims = 1:num_pcs, resolution = resolution);
108 } else {
109 seuset <- Seurat::RunTSNE(seuset, dims = 1:num_pcs, resolution = resolution, perplexity = perplexity);
110 }
102 Seurat::DimPlot(seuset, reduction = "tsne") 111 Seurat::DimPlot(seuset, reduction = "tsne")
103 112
104 #+ echo = FALSE 113 #+ echo = FALSE
105 if(showcode == TRUE && heatmaps == TRUE) print("Marker Genes") 114 if (showcode == TRUE && heatmaps == TRUE) print("Marker Genes")
106 #+ echo = `showcode`, warning = `warn`, include = `heatmaps` 115 #+ echo = `showcode`, warning = `warn`, include = `heatmaps`
107 markers <- Seurat::FindAllMarkers(seuset, only.pos = TRUE, min.pct = min_pct, logfc.threshold = logfc_threshold) 116 markers <- Seurat::FindAllMarkers(seuset, only.pos = TRUE, min.pct = min_pct, logfc.threshold = logfc_threshold)
108 top10 <- dplyr::group_by(markers, cluster) 117 top10 <- dplyr::group_by(markers, cluster)
109 top10 <- dplyr::top_n(top10, 10, avg_logFC) 118 top10 <- dplyr::top_n(top10, 10, avg_logFC)
110 Seurat::DoHeatmap(seuset, features = top10$gene) 119 Seurat::DoHeatmap(seuset, features = top10$gene)
120 # nolint end