9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
|
# File 'lib/macroape/cli/preprocess_collection.rb', line 9
def self.main(argv)
doc = " Command-line format:\n \#{run_tool_cmd} <file or folder with PWMs or .stdin with filenames> <output file> [options]\n\n Options:\n [-p <list of P-values>] - comma separated(no spaces allowed) list of P-values to precalculate thresholds\n [-d <rough discretization>,<precise discretization>] - set discretization rates, comma delimited (no spaces allowed), order doesn't matter\n [--silent] - hide current progress information during scan (printed to stderr by default)\n [--pcm] - treat the input file as Position Count Matrix. PCM-to-PWM transformation to be done internally.\n [--boundary lower|upper] Upper boundary (default) means that the obtained P-value is greater than or equal to the requested P-value\n [-b <background probabilities] ACGT - 4 numbers, comma-delimited(spaces not allowed), sum should be equal to 1, like 0.25,0.24,0.26,0.25\n\n The tool preprocesses and stores Macroape motif collection in the specified YAML-file.\n\n Example:\n \#{run_tool_cmd} ./motifs collection.yaml -p 0.001,0.0005,0.0001 -d 1,10 -b 0.2,0.3,0.3,0.2\n EOS\n\n if argv.empty? || ['-h', '--h', '-help', '--help'].any?{|help_option| argv.include?(help_option)}\n $stderr.puts doc\n exit\n end\n\n data_model = argv.delete('--pcm') ? Bioinform::PCM : Bioinform::PWM\n\n default_pvalues = [0.0005]\n background = [1,1,1,1]\n rough_discretization = 1\n precise_discretization = 10\n max_hash_size = 10000000\n\n data_source = argv.shift\n output_file = argv.shift\n\n raise 'No input. You should specify file or folder with pwms' unless data_source\n raise \"Error! File or folder \#{data_source} doesn't exist\" unless Dir.exist?(data_source) || File.exist?(data_source) || data_source == '.stdin'\n raise 'You should specify output file' unless output_file\n\n pvalues = []\n silent = false\n pvalue_boundary = :upper\n\n until argv.empty?\n case argv.shift\n when '-b'\n background = argv.shift.split(',').map(&:to_f)\n raise 'background should be symmetric: p(A)=p(T) and p(G) = p(C)' unless background == background.reverse\n when '-p'\n pvalues = argv.shift.split(',').map(&:to_f)\n when '-d'\n rough_discretization, precise_discretization = argv.shift.split(',').map(&:to_f).sort\n when '--max-hash-size'\n max_hash_size = argv.shift.to_i\n when '--silent'\n silent = true\n when '--boundary'\n pvalue_boundary = argv.shift.to_sym\n raise 'boundary should be either lower or upper' unless pvalue_boundary == :lower || pvalue_boundary == :upper\n end\n end\n pvalues = default_pvalues if pvalues.empty?\n\n collection = Bioinform::Collection.new(rough_discretization: rough_discretization,\n precise_discretization: precise_discretization,\n background: background,\n pvalues: pvalues)\n\n data_source = data_source.gsub(\"\\\\\",'/')\n if File.directory?(data_source)\n motifs = Dir.glob(File.join(data_source,'*')).sort.map do |filename|\n pwm = data_model.new(File.read(filename))\n pwm.name ||= File.basename(filename, File.extname(filename))\n pwm\n end\n elsif File.file?(data_source)\n input = File.read(data_source)\n motifs = data_model.split_on_motifs(input)\n elsif data_source == '.stdin'\n filelist = $stdin.read.shellsplit\n motifs = []\n filelist.each do |filename|\n motif = data_model.new(File.read(filename))\n motif.name ||= File.basename(filename, File.extname(filename))\n motif.set_parameters(background: background)\n motifs << motif\n end\n else\n raise \"Specified data source `\#{data_source}` is neither directory nor file nor even .stdin\"\n end\n\n pwms = motifs.map(&:to_pwm)\n\n pwms.each_with_index do |pwm,index|\n $stderr.puts \"Motif \#{pwm.name}, length: \#{pwm.length} (\#{index+1} of \#{pwms.size}, \#{index*100/pwms.size}% complete)\" unless silent\n\n # When support of onefile collections is introduced - then here should be check if name exists.\n # Otherwise it should skip motif and tell you about this\n # Also two command line options to fail on skipping or to skip silently should be included\n\n info = OpenStruct.new(rough: {}, precise: {})\n pwm.set_parameters(background: background, max_hash_size: max_hash_size)\n skip_motif = false\n\n\n fill_rough_infos = ->(pvalue, threshold, real_pvalue) do\n if real_pvalue == 0\n $stderr.puts \"\#{pwm.name} at pvalue \#{pvalue} has threshold that yields real-pvalue 0 in rough mode. Rough calculation will be skipped\"\n else\n info.rough[pvalue] = threshold / rough_discretization\n end\n end\n\n fill_precise_infos = ->(pvalue, threshold, real_pvalue) do\n if real_pvalue == 0\n $stderr.puts \"\#{pwm.name} at pvalue \#{pvalue} has threshold that yields real-pvalue 0 in precise mode. Motif will be excluded from collection\"\n skip_motif = true\n else\n info.precise[pvalue] = threshold / precise_discretization\n end\n end\n\n if pvalue_boundary == :lower\n pwm.discrete(rough_discretization).thresholds(*pvalues, &fill_rough_infos)\n else\n pwm.discrete(rough_discretization).weak_thresholds(*pvalues, &fill_rough_infos)\n end\n\n if pvalue_boundary == :lower\n pwm.discrete(precise_discretization).thresholds(*pvalues, &fill_precise_infos)\n else\n pwm.discrete(precise_discretization).weak_thresholds(*pvalues,&fill_precise_infos)\n end\n collection.add_pm(pwm, info) unless skip_motif\n end\n $stderr.puts \"100% complete. Saving results\" unless silent\n File.open(output_file, 'w') do |f|\n f.puts(collection.to_yaml)\n end\n puts OutputInformation.new{|infos|\n infos.add_parameter('P', 'P-value list', pvalues.join(','))\n infos.add_parameter('VR', 'discretization value, rough', rough_discretization)\n infos.add_parameter('VP', 'discretization value, precise', precise_discretization)\n infos.add_parameter('PB', 'P-value boundary', pvalue_boundary)\n infos.background_parameter('B', 'background', background)\n }.result\nrescue => err\n $stderr.puts \"\\n\#{err}\\n\#{err.backtrace.first(5).join(\"\\n\")}\\n\\nUse --help option for help\\n\\n\#{doc}\"\nend\n".strip_doc
|