|
0
|
1 ###############################################################################
|
|
|
2 # title: chipenrich
|
|
|
3 # author: Xiaowei
|
|
|
4 # time: Mar.31 2021
|
|
|
5 ###############################################################################
|
|
|
6
|
|
|
7 spec <- matrix(c("input_peaks", "P",1, "character", "peaks file",
|
|
|
8 "input_genome","G",1, "character", "genome",
|
|
|
9 "input_locusdef","L",1,"character", "locusdef",
|
|
|
10 "input_mappability","Map",0, "logical", "mappability",
|
|
|
11 "input_qc_plot","Q",0, "logical", "qc_plot",
|
|
|
12 "input_method","M",1,"character", "method",
|
|
|
13 "input_geneset","GS",1,"character", "geneset_category",
|
|
|
14 "input_minSize","Min",1,"integer","minSize",
|
|
|
15 "input_maxSize","Max",1,"integer", "maxSize",
|
|
|
16 "input_randomization","random",1,"character","randomization",
|
|
|
17 "input_reglocation","reg",1,"character","reglocation",
|
|
|
18 "input_num_peak_threshould","threshould",1,"numeric", "peak_threshould",
|
|
|
19 "output_peaks","out_peaks",1,"character","output_peaks",
|
|
|
20 "output_peaks_per_gene","out_ppg",1,"character","peaks per gene",
|
|
|
21 "output_enrich_results","out_er",1,"character", "enrich result",
|
|
|
22 "output_qc_plot", "out_qc", 1,"character","output_qc_plot"
|
|
|
23 ),byrow = TRUE, ncol = 5)
|
|
|
24
|
|
|
25 opt <- getopt::getopt(spec)
|
|
|
26
|
|
|
27
|
|
|
28 #----------------------------------------------------------------
|
|
|
29
|
|
|
30 input_peaks <- read.csv(opt$input_peaks,header = TRUE)
|
|
|
31 input_genome = opt$input_genome
|
|
|
32 input_locusdef = opt$input_locusdef
|
|
|
33
|
|
|
34 if (!is.null(opt$input_mappability)){
|
|
|
35 input_mappability = 24 #根据表格来填写
|
|
|
36 }else{
|
|
|
37 input_mappability = NULL #supported_read_lengths()
|
|
|
38 }
|
|
|
39
|
|
|
40 if(!is.null(opt$input_qc_plot)){
|
|
|
41 input_qc_plot = opt$input_qc_plot
|
|
|
42 input_output_name = "x"
|
|
|
43 }else{
|
|
|
44 input_qc_plot = FALSE
|
|
|
45 input_output_name = NULL
|
|
|
46 }
|
|
|
47
|
|
|
48
|
|
|
49 if (!is.null(opt$input_method)){input_method = opt$input_method}else{input_method ="chipenrich"}
|
|
|
50 if(!is.null(opt$input_geneset)){input_geneset = strsplit(opt$input_geneset, split = ",")[[1]]}else{input_geneset = NULL}
|
|
|
51 if(!is.null(opt$input_minSize)){input_minSize = opt$input_minSize}else{input_minSize = 15}
|
|
|
52 if(!is.null(opt$input_maxSize)){input_maxSize = opt$input_maxSize}else{input_maxSize = 2000}
|
|
|
53 if(!is.null(opt$input_randomization) & opt$input_randomization != "NULL"){input_randomization = opt$input_randomization}else{input_randomization = NULL}
|
|
|
54 input_core = parallel::detectCores()
|
|
|
55 if(!is.null(opt$input_reglocation)){input_reglocation = opt$input_reglocation}else{input_reglocation = "tss"} #proxReg
|
|
|
56 #chipenrich、hybridenrich
|
|
|
57 if(!is.null(opt$input_num_peak_threshould)){input_num_peak_threshould = opt$input_num_peak_threshould}else{input_num_peak_threshould = 1}
|
|
|
58 #================================================================
|
|
|
59 #run codes
|
|
|
60 #================================================================
|
|
|
61 suppressPackageStartupMessages(library(chipenrich))
|
|
|
62 ## ---- chipenrich----------------------------------------
|
|
|
63 if (input_method == "chipenrich"){
|
|
|
64 results = chipenrich(peaks = input_peaks, genome = input_genome, genesets = input_geneset,
|
|
|
65 locusdef = input_locusdef, qc_plots = FALSE, out_name = input_output_name,
|
|
|
66 mappability = input_mappability,
|
|
|
67 min_geneset_size = input_minSize,
|
|
|
68 max_geneset_size = input_maxSize,
|
|
|
69 randomization = input_randomization,
|
|
|
70 num_peak_threshold = input_num_peak_threshould,
|
|
|
71 n_cores = input_core)
|
|
|
72 }
|
|
|
73
|
|
|
74
|
|
|
75 ## ----polyenrich----------------------------------------
|
|
|
76 if (input_method == "polyenrich"){
|
|
|
77 results = polyenrich(peaks = input_peaks, genome = input_genome, genesets = input_geneset,
|
|
|
78 method = 'polyenrich',
|
|
|
79 locusdef = input_locusdef, qc_plots = FALSE, out_name = input_output_name,
|
|
|
80 mappability = input_mappability,
|
|
|
81 min_geneset_size = input_minSize,
|
|
|
82 max_geneset_size = input_maxSize,
|
|
|
83 randomization = input_randomization,
|
|
|
84 n_cores = input_core)
|
|
|
85 }
|
|
|
86
|
|
|
87 ## ----hybridenrich----------------------------------------
|
|
|
88 if (input_method == "hybridenrich"){
|
|
|
89 results = hybridenrich(peaks = input_peaks, genome = input_genome, genesets = input_geneset,
|
|
|
90 locusdef = input_locusdef, qc_plots = FALSE, out_name = input_output_name,
|
|
|
91 mappability = input_mappability,
|
|
|
92 min_geneset_size = input_minSize,
|
|
|
93 max_geneset_size = input_maxSize,
|
|
|
94 randomization = input_randomization,
|
|
|
95 num_peak_threshold = input_num_peak_threshould,
|
|
|
96 n_cores = input_core)
|
|
|
97 }
|
|
|
98
|
|
|
99 ## ----proxReg----------------------------------------
|
|
|
100 if (input_method == "proxReg"){
|
|
|
101 results = proxReg(peaks = input_peaks, reglocation = input_reglocation,
|
|
|
102 genome = input_genome, genesets=input_geneset,
|
|
|
103 min_geneset_size = input_minSize,
|
|
|
104 max_geneset_size = input_maxSize,
|
|
|
105 randomization = input_randomization,
|
|
|
106 out_name=input_output_name,
|
|
|
107 n_cores = input_core)
|
|
|
108 }
|
|
|
109
|
|
|
110 ## ----broadenrich----------------------------------------peaks_H3K4me3_GM12878---
|
|
|
111 if (input_method == "broadenrich"){
|
|
|
112 results = broadenrich(peaks = input_peaks, genome = input_genome, genesets = input_geneset,
|
|
|
113 locusdef = input_locusdef, qc_plots = FALSE, out_name = input_output_name,
|
|
|
114 mappability = input_mappability,
|
|
|
115 min_geneset_size = input_minSize,
|
|
|
116 max_geneset_size = input_maxSize,
|
|
|
117 randomization = input_randomization,
|
|
|
118 n_cores = input_core)
|
|
|
119 }
|
|
|
120
|
|
|
121 #===============================================================================
|
|
|
122 # 输出
|
|
|
123 #===============================================================================
|
|
|
124 output_peaks = opt$output_peaks
|
|
|
125 output_peaks_per_gene = opt$output_peaks_per_gene
|
|
|
126 output_enrich_result = opt$output_enrich_results
|
|
|
127
|
|
|
128 if (!is.null(output_peaks)){
|
|
|
129 write.csv(results$peaks, file = output_peaks)
|
|
|
130 }
|
|
|
131
|
|
|
132 if (!is.null(output_peaks_per_gene)){
|
|
|
133 if (input_method != "proxReg"){
|
|
|
134 write.csv(results$peaks_per_gene, file = output_peaks_per_gene)
|
|
|
135 }
|
|
|
136 }
|
|
|
137
|
|
|
138 if (!is.null(output_enrich_result)){
|
|
|
139 write.csv(results$results, file = output_enrich_result)
|
|
|
140 }
|
|
|
141
|
