Mercurial > repos > gaelcge > r_signac_galaxy

comparison signac-runQC.R @ 0:6e0b320d8b6a draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit dc808171975d0012e25bd7b32adc7a5a5c56a145-dirty"
author gaelcge
date Tue, 02 Aug 2022 19:11:27 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6e0b320d8b6a
1 #!/usr/bin/env Rscript
2
3 # Load optparse we need to check inputs
4
5 suppressPackageStartupMessages(require(optparse))
6
7 # Load common functions
8
9 suppressPackageStartupMessages(require(workflowscriptscommon))
10
11 # parse options
12
13 option_list = list(
14 make_option(
15 c("--signac-object"),
16 action = "store",
17 default = NA,
18 type = 'character',
19 help = ""
20 ),
21 make_option(
22 c("--tss-threshold"),
23 action = "store",
24 default = NA,
25 type = 'character',
26 help = "TSS enrichment threshold for marking regions as high tss regions."
27 ),
28 make_option(
29 c("--output-tss-plot"),
30 action = "store",
31 default = NA,
32 type = 'character',
33 help = "TSS output plot."
34 ),
35 make_option(
36 c("--frag-history-plot"),
37 action = "store",
38 default = NA,
39 type = 'character',
40 help = "Fragment length periodicity plot."
41 ),
42 make_option(
43 c("--fragment-file"),
44 action = "store",
45 default = NA,
46 type = 'character',
47 help = "Fragment file."
48 ),
49 make_option(
50 c("-w", "--png-width"),
51 action = "store",
52 default = 1000,
53 type = 'integer',
54 help = "Width of png (px)."
55 ),
56 make_option(
57 c("-j", "--png-height"),
58 action = "store",
59 default = 1000,
60 type = 'integer',
61 help = "Height of png file (px)."
62 ),
63 make_option(
64 c("--output-object-file"),
65 action = "store",
66 default = NA,
67 type = 'character',
68 help = "File name in which to store serialized R matrix object."
69 )
70 )
71
72 opt <- wsc_parse_args(option_list)
73
74 suppressPackageStartupMessages(require(Seurat))
75 suppressPackageStartupMessages(require(Signac))
76 suppressPackageStartupMessages(require(GenomeInfoDb))
77 suppressPackageStartupMessages(require(EnsDb.Hsapiens.v75))
78 suppressPackageStartupMessages(require(EnsDb.Mmusculus.v79))
79
80 set.seed(1234)
81
82 # extract gene annotations from EnsDb
83 signac_object <- readRDS(file = opt$signac_object)
84
85 ## Modify fragments file location
86 # current_wd <- getwd()
87 # new_framgent_file_loc <- paste(current_wd,"fragments.tsv.gz",sep = "/")
88 # signac_object@assays$peaks@fragments[[1]]@path <- new_framgent_file_loc
89 # signac_object@assays$peaks@fragments[[1]]@path <- opt$fragment_file
90 #print(normalizePath(path = paste0(opt$fragment_file, '.tbi'), mustWork = TRUE))
91 #system(paste0("mv ", paste0(opt$fragment_file, '.tbi'), ' ', normalizePath(path = paste0(opt$fragment_file, '.tbi'), mustWork = TRUE)))
92 # The error was that the fragment file and fragment index file were not being imported into galaxy properly. They should be imported as tabluar.gz
93 signac_object@assays$peaks@fragments[[1]] <- UpdatePath(signac_object@assays$peaks@fragments[[1]], opt$fragment_file)
94
95 # compute nucleosome signal score per cell
96 signac_object <- NucleosomeSignal(object = signac_object)
97
98 # compute TSS enrichment score per cell
99 signac_object <- TSSEnrichment(object = signac_object, fast = FALSE)
100
101 # add blacklist ratio and fraction of reads in peaks
102 signac_object$pct_reads_in_peaks <- signac_object$peak_region_fragments / signac_object$passed_filters * 100
103 signac_object$blacklist_ratio <- signac_object$blacklist_region_fragments / signac_object$peak_region_fragments
104
105 signac_object$high.tss <- ifelse(signac_object$TSS.enrichment > as.numeric(opt$tss_threshold), 'High', 'Low')
106
107 png(filename = opt$output_tss_plot, width = opt$png_width, height = opt$png_height)
108 TSSPlot(signac_object, group.by = 'high.tss') + NoLegend()
109 dev.off()
110
111 signac_object$nucleosome_group <- ifelse(signac_object$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
112
113 png(filename = opt$frag_history_plot, width = opt$png_width, height = opt$png_height)
114 FragmentHistogram(object = signac_object, group.by = 'nucleosome_group')
115 dev.off()
116
117 print(summary(signac_object@meta.data))
118
119 # Output to a serialized R object
120 saveRDS(signac_object, file = opt$output_object_file)