Mercurial > repos > gaelcge > r_signac_galaxy
comparison signac-runQC.R @ 0:6e0b320d8b6a draft default tip
"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) |