Mercurial > repos > azomics > metacyto_histogram

comparison metacyto_histogram.R @ 0:f5526d97056c draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/AstraZeneca-Omics/immport-galaxy-tools/tree/master/flowtools/metacyto_histogram commit cb978232e32b64f7b0ff3c1852e708361045d268"
author azomics
date Thu, 29 Jul 2021 22:15:11 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f5526d97056c
1 #!/usr/bin/env Rscript
2 ######################################################################
3 # Copyright (c) 2018 Northrop Grumman.
4 # All rights reserved.
5 ######################################################################
6 #
7 # Version 1 - January 2018
8 # Author: Cristel Thomas
9 #
10 #
11
12 library(flowCore)
13 library(MetaCyto)
14
15 check_cluster_def <- function(cl_def) {
16 if (cl_def == "" || cl_def == "None") {
17 quit(save = "no", status = 12, runLast = FALSE)
18 } else {
19 tmp <- gsub(" ", "", cl_def, fixed = TRUE)
20 clean_def <- gsub(",", "|", tmp, fixed = TRUE)
21 return(clean_def)
22 }
23 }
24
25 generate_plots <- function(fpath = "", fname = "", gates = vector(), outdir = "", uc = "",
26 flag_pdf = F) {
27 dir.create(outdir)
28 ff <- read.FCS(fpath, truncate_max_range = F)
29 markers <- markerFinder(ff)
30 colnames(ff@exprs) <- markers
31
32 sc <- searchCluster(fcsFrame = ff, clusterLabel = gates)
33
34 if (length(gates) == length(sc$clusterList)) {
35 sink(uc)
36 cat("All provided cluster definition were used.")
37 sink()
38 } else {
39 unused_cluster <- setdiff(gates, names(sc$clusterList))
40 write.table(unused_cluster, uc, quote = F, row.names = F, col.names = F)
41 }
42
43 groupname <- unlist(strsplit(fname, ".fcs"))[[1]]
44 extension <- if (flag_pdf) "plot.pdf" else "plot.png"
45 for (i in seq_len(length(sc$clusterList))) {
46 gate <- gsub("|", "", names(sc$clusterList[i]), fixed = T)
47 plotname <- paste(c(groupname, gate, extension), collapse = "_")
48 outplot <- file.path(outdir, plotname)
49 if (flag_pdf) {
50 pdf(outplot, useDingbats = F, onefile = T)
51 par(mfrow = c(2, 2))
52 for (j in seq_len(length(markers))) {
53 if (markers[[j]] != "SAMPLE_ID" && markers[[j]] != "TIME") {
54 plot_title <- paste0(markers[[j]], ", cluster definition:\n", gate)
55 x_all <- ff@exprs[, markers[[j]]]
56 b <- seq(min(x_all), max(x_all), ((max(x_all) - min(x_all)) / 100))
57 subset <- ff@exprs[sc$clusterList[[i]], markers[[j]]]
58 hist(x_all, col = rgb(0, 0, 0), xlab = markers[[j]], breaks = b, freq = T,
59 border = F, main = plot_title)
60 hist(subset, add = T, breaks = b, col = rgb(1, 0, 0), freq = T, border = F)
61 if (markers[[j]] %in% names(sc$cutoff)) {
62 abline(v = sc$cutoff[markers[[j]]])
63 }
64 }
65 }
66 dev.off()
67 } else {
68 markers_ct <- length(markers) - length(grep(x = markers, pattern = "SAMPLE_ID|TIME"))
69 nb_rows <- ceiling(markers_ct / 2)
70 h <- nb_rows * 400
71 png(outplot, type = "cairo", height = h, width = 800)
72 par(mfrow = c(nb_rows, 2))
73 for (j in seq_len(length(markers))) {
74 if (markers[[j]] != "SAMPLE_ID" && markers[[j]] != "TIME") {
75 plot_title <- paste0(markers[[j]], ", cluster definition:\n", gate)
76 x_all <- ff@exprs[, markers[[j]]]
77 b <- seq(min(x_all), max(x_all), ((max(x_all) - min(x_all)) / 100))
78 subset <- ff@exprs[sc$clusterList[[i]], markers[[j]]]
79 hist(x_all, col = rgb(0, 0, 0), xlab = markers[[j]], breaks = b, freq = T,
80 border = F, main = plot_title)
81 hist(subset, add = T, breaks = b, col = rgb(1, 0, 0), freq = T, border = F)
82 if (markers[[j]] %in% names(sc$cutoff)) {
83 abline(v = sc$cutoff[markers[[j]]])
84 }
85 }
86 }
87 dev.off()
88 }
89 }
90 }
91
92 check_fcs_file <- function(inputf="", inputn="", clusters=vector(),
93 output_dir = "", unused = "", flag = F) {
94 is_valid <- FALSE
95 tryCatch({
96 is_valid <- isFCSfile(inputf)
97 }, error = function(ex) {
98 print(paste("Input file is not a valid FCS file.", ex))
99 })
100 if (is_valid) {
101 generate_plots(inputf, inputn, clusters, output_dir, unused, flag)
102 } else {
103 quit(save = "no", status = 12, runLast = FALSE)
104 }
105 }
106
107 ################################################################################
108 ################################################################################
109 args <- commandArgs(trailingOnly = TRUE)
110
111 gates <- vector()
112 if (args[6] == "F") {
113 ## obvs deal with it if file
114 cluster_file <- read.table(args[7], header = F, colClasses = "character")
115 gates <- unlist(cluster_file)
116 } else {
117 cl_df <- args[7:length(args)]
118 gates <- sapply(cl_df, check_cluster_def)
119 }
120
121 flag_pdf <- if (args[5] == "PDF") TRUE else FALSE
122 gate_list <- toupper(gates)
123 check_fcs_file(args[1], args[2], gate_list, args[3], args[4], flag_pdf)