Mercurial > repos > rnateam > graphprot_predict_profile
comparison graphprot_predict_wrapper.py @ 5:58ebf089377e draft default tip
planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/tools/rna_tools/graphprot commit 902e994cb04db968ce797afa676f5fa6512ab6d3
| author | bgruening |
|---|---|
| date | Tue, 06 Aug 2024 14:54:58 +0000 |
| parents | 2d6de5b769e6 |
| children |
comparison
equal
deleted
inserted
replaced
| 4:2d6de5b769e6 | 5:58ebf089377e |
|---|---|
| 4 import os | 4 import os |
| 5 import subprocess | 5 import subprocess |
| 6 import sys | 6 import sys |
| 7 | 7 |
| 8 import gplib | 8 import gplib |
| 9 | |
| 10 | 9 |
| 11 """ | 10 """ |
| 12 | 11 |
| 13 TOOL DEPENDENCIES | 12 TOOL DEPENDENCIES |
| 14 ================= | 13 ================= |
| 79 --sc-thr 0.0 --max-merge-dist 0 --conf-out --ap-extlr 5 | 78 --sc-thr 0.0 --max-merge-dist 0 --conf-out --ap-extlr 5 |
| 80 | 79 |
| 81 """ | 80 """ |
| 82 | 81 |
| 83 | 82 |
| 84 ############################################################################### | 83 ####################################################################### |
| 84 | |
| 85 | 85 |
| 86 def setup_argument_parser(): | 86 def setup_argument_parser(): |
| 87 """Setup argparse parser.""" | 87 """Setup argparse parser.""" |
| 88 help_description = """ | 88 help_description = """ |
| 89 Galaxy wrapper script for GraphProt (-action predict and -action | 89 Galaxy wrapper script for GraphProt (-action predict and -action |
| 98 If --gen-site-bed .bed file is provided, peak regions will be output | 98 If --gen-site-bed .bed file is provided, peak regions will be output |
| 99 with genomic coordinates too. | 99 with genomic coordinates too. |
| 100 | 100 |
| 101 """ | 101 """ |
| 102 # Define argument parser. | 102 # Define argument parser. |
| 103 p = ap.ArgumentParser(add_help=False, | 103 p = ap.ArgumentParser( |
| 104 prog="graphprot_predict_wrapper.py", | 104 add_help=False, |
| 105 description=help_description, | 105 prog="graphprot_predict_wrapper.py", |
| 106 formatter_class=ap.MetavarTypeHelpFormatter) | 106 description=help_description, |
| 107 formatter_class=ap.MetavarTypeHelpFormatter, | |
| 108 ) | |
| 107 | 109 |
| 108 # Argument groups. | 110 # Argument groups. |
| 109 p_man = p.add_argument_group("REQUIRED ARGUMENTS") | 111 p_man = p.add_argument_group("REQUIRED ARGUMENTS") |
| 110 p_opt = p.add_argument_group("OPTIONAL ARGUMENTS") | 112 p_opt = p.add_argument_group("OPTIONAL ARGUMENTS") |
| 111 | 113 |
| 112 # Required arguments. | 114 # Required arguments. |
| 113 p_opt.add_argument("-h", "--help", | 115 p_opt.add_argument("-h", "--help", action="help", help="Print help message") |
| 114 action="help", | 116 p_man.add_argument( |
| 115 help="Print help message") | 117 "--fasta", |
| 116 p_man.add_argument("--fasta", | 118 dest="in_fa", |
| 117 dest="in_fa", | 119 type=str, |
| 118 type=str, | 120 required=True, |
| 119 required=True, | 121 help="Sequences .fa file to predict" " on (option -fasta)", |
| 120 help="Sequences .fa file to predict" | 122 ) |
| 121 " on (option -fasta)") | 123 p_man.add_argument( |
| 122 p_man.add_argument("--model", | 124 "--model", |
| 123 dest="in_model", | 125 dest="in_model", |
| 124 type=str, | 126 type=str, |
| 125 required=True, | 127 required=True, |
| 126 help="GraphProt model file to use for predictions" | 128 help="GraphProt model file to use for predictions" " (option -model)", |
| 127 " (option -model)") | 129 ) |
| 128 p_man.add_argument("--params", | 130 p_man.add_argument( |
| 129 dest="in_params", | 131 "--params", |
| 130 type=str, | 132 dest="in_params", |
| 131 required=True, | 133 type=str, |
| 132 help="Parameter file for given model") | 134 required=True, |
| 133 p_man.add_argument("--data-id", | 135 help="Parameter file for given model", |
| 134 dest="data_id", | 136 ) |
| 135 type=str, | 137 p_man.add_argument( |
| 136 required=True, | 138 "--data-id", |
| 137 help="Data ID (option -prefix)") | 139 dest="data_id", |
| 140 type=str, | |
| 141 required=True, | |
| 142 help="Data ID (option -prefix)", | |
| 143 ) | |
| 138 # ---> I'm a conditional argument <--- | 144 # ---> I'm a conditional argument <--- |
| 139 p_opt.add_argument("--ws-pred", | 145 p_opt.add_argument( |
| 140 dest="ws_pred", | 146 "--ws-pred", |
| 141 default=False, | 147 dest="ws_pred", |
| 142 action="store_true", | 148 default=False, |
| 143 help="Run a whole site prediction instead " | 149 action="store_true", |
| 144 "of calculating profiles (default: false)") | 150 help="Run a whole site prediction instead " |
| 151 "of calculating profiles (default: false)", | |
| 152 ) | |
| 145 # Additional arguments. | 153 # Additional arguments. |
| 146 p_opt.add_argument("--sc-thr", | 154 p_opt.add_argument( |
| 147 dest="score_thr", | 155 "--sc-thr", |
| 148 type=float, | 156 dest="score_thr", |
| 149 default=0, | 157 type=float, |
| 150 help="Score threshold for extracting " | 158 default=0, |
| 151 "average profile peak regions (default: 0)") | 159 help="Score threshold for extracting " |
| 152 p_opt.add_argument("--max-merge-dist", | 160 "average profile peak regions (default: 0)", |
| 153 dest="max_merge_dist", | 161 ) |
| 154 type=int, | 162 p_opt.add_argument( |
| 155 default=0, | 163 "--max-merge-dist", |
| 156 choices=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], | 164 dest="max_merge_dist", |
| 157 help="Maximum merge distance for nearby peak regions" | 165 type=int, |
| 158 " (default: report all non-overlapping regions)") | 166 default=0, |
| 159 p_opt.add_argument("--gen-site-bed", | 167 choices=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], |
| 160 dest="genomic_sites_bed", | 168 help="Maximum merge distance for nearby peak regions" |
| 161 type=str, | 169 " (default: report all non-overlapping regions)", |
| 162 help=".bed file specifying the genomic regions of " | 170 ) |
| 163 "the input .fa sequences. Corrupt .bed " | 171 p_opt.add_argument( |
| 164 "information will be punished (default: false)") | 172 "--gen-site-bed", |
| 165 p_opt.add_argument("--conf-out", | 173 dest="genomic_sites_bed", |
| 166 dest="conf_out", | 174 type=str, |
| 167 default=False, | 175 help=".bed file specifying the genomic regions of " |
| 168 action="store_true", | 176 "the input .fa sequences. Corrupt .bed " |
| 169 help="Output filtered peak regions BED file or " | 177 "information will be punished (default: false)", |
| 170 "predictions file (if --ws-pred) using the median " | 178 ) |
| 171 "positive training site score for filtering " | 179 p_opt.add_argument( |
| 172 "(default: false)") | 180 "--conf-out", |
| 173 p_opt.add_argument("--gp-output", | 181 dest="conf_out", |
| 174 dest="gp_output", | 182 default=False, |
| 175 default=False, | 183 action="store_true", |
| 176 action="store_true", | 184 help="Output filtered peak regions BED file or " |
| 177 help="Print output produced by GraphProt " | 185 "predictions file (if --ws-pred) using the median " |
| 178 "(default: false)") | 186 "positive training site score for filtering " |
| 179 p_opt.add_argument("--ap-extlr", | 187 "(default: false)", |
| 180 dest="ap_extlr", | 188 ) |
| 181 type=int, | 189 p_opt.add_argument( |
| 182 default=5, | 190 "--gp-output", |
| 183 choices=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], | 191 dest="gp_output", |
| 184 help="Define average profile up- and downstream " | 192 default=False, |
| 185 "extension to produce the average profile. " | 193 action="store_true", |
| 186 "The mean over small sequence windows " | 194 help="Print output produced by GraphProt " "(default: false)", |
| 187 "(window length = --ap-extlr*2 + 1) is used to " | 195 ) |
| 188 "get position scores, thus the average profile " | 196 p_opt.add_argument( |
| 189 "is more smooth than the initial profile output " | 197 "--ap-extlr", |
| 190 "by GraphProt (default: 5)") | 198 dest="ap_extlr", |
| 199 type=int, | |
| 200 default=5, | |
| 201 choices=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], | |
| 202 help="Define average profile up- and downstream " | |
| 203 "extension to produce the average profile. " | |
| 204 "The mean over small sequence windows " | |
| 205 "(window length = --ap-extlr*2 + 1) is used to " | |
| 206 "get position scores, thus the average profile " | |
| 207 "is more smooth than the initial profile output " | |
| 208 "by GraphProt (default: 5)", | |
| 209 ) | |
| 191 return p | 210 return p |
| 192 | 211 |
| 193 | 212 |
| 194 ############################################################################### | 213 ####################################################################### |
| 195 | 214 |
| 196 if __name__ == '__main__': | 215 if __name__ == "__main__": |
| 197 | 216 |
| 198 # Setup argparse. | 217 # Setup argparse. |
| 199 parser = setup_argument_parser() | 218 parser = setup_argument_parser() |
| 200 # Read in command line arguments. | 219 # Read in command line arguments. |
| 201 args = parser.parse_args() | 220 args = parser.parse_args() |
| 207 # Check for Linux. | 226 # Check for Linux. |
| 208 assert "linux" in sys.platform, "please use Linux" | 227 assert "linux" in sys.platform, "please use Linux" |
| 209 # Check tool availability. | 228 # Check tool availability. |
| 210 assert gplib.is_tool("GraphProt.pl"), "GraphProt.pl not in PATH" | 229 assert gplib.is_tool("GraphProt.pl"), "GraphProt.pl not in PATH" |
| 211 # Check file inputs. | 230 # Check file inputs. |
| 212 assert os.path.exists(args.in_fa), \ | 231 assert os.path.exists(args.in_fa), 'input .fa file "%s" not found' % (args.in_fa) |
| 213 "input .fa file \"%s\" not found" % (args.in_fa) | 232 assert os.path.exists(args.in_model), 'input .model file "%s" not found' % ( |
| 214 assert os.path.exists(args.in_model), \ | 233 args.in_model |
| 215 "input .model file \"%s\" not found" % (args.in_model) | 234 ) |
| 216 assert os.path.exists(args.in_params), \ | 235 assert os.path.exists(args.in_params), 'input .params file "%s" not found' % ( |
| 217 "input .params file \"%s\" not found" % (args.in_params) | 236 args.in_params |
| 237 ) | |
| 218 # Count .fa entries. | 238 # Count .fa entries. |
| 219 c_in_fa = gplib.count_fasta_headers(args.in_fa) | 239 c_in_fa = gplib.count_fasta_headers(args.in_fa) |
| 220 assert c_in_fa, "input .fa file \"%s\" no headers found" % (args.in_fa) | 240 assert c_in_fa, 'input .fa file "%s" no headers found' % (args.in_fa) |
| 221 print("# input .fa sequences: %i" % (c_in_fa)) | 241 print("# input .fa sequences: %i" % (c_in_fa)) |
| 222 # Read in FASTA sequences to check for uppercase sequences. | 242 # Read in FASTA sequences to check for uppercase sequences. |
| 223 seqs_dic = gplib.read_fasta_into_dic(args.in_fa) | 243 seqs_dic = gplib.read_fasta_into_dic(args.in_fa) |
| 224 # Check for lowercase only sequences, which cause GP to crash. | 244 # Check for lowercase only sequences, which cause GP to crash. |
| 225 error_mess = "input sequences encountered containing "\ | 245 error_mess = ( |
| 226 "only lowercase characters or lowercase characters in between "\ | 246 "input sequences encountered containing " |
| 227 "uppercase characters. Please provide either all uppercase "\ | 247 "only lowercase characters or lowercase characters in between " |
| 228 "sequences or sequences containing uppercase regions surrounded "\ | 248 "uppercase characters. Please provide either all uppercase " |
| 229 "by lowercase context regions for structure calculation (see "\ | 249 "sequences or sequences containing uppercase regions surrounded " |
| 230 "viewpoint concept in original GraphProt publication "\ | 250 "by lowercase context regions for structure calculation (see " |
| 251 "viewpoint concept in original GraphProt publication " | |
| 231 "for more details)" | 252 "for more details)" |
| 253 ) | |
| 232 if args.ws_pred: | 254 if args.ws_pred: |
| 233 bad_ids = gplib.check_seqs_dic_format(seqs_dic) | 255 bad_ids = gplib.check_seqs_dic_format(seqs_dic) |
| 234 assert not bad_ids, "%s" % (error_mess) | 256 assert not bad_ids, "%s" % (error_mess) |
| 235 | 257 |
| 236 c_uc_nt = gplib.seqs_dic_count_uc_nts(seqs_dic) | 258 c_uc_nt = gplib.seqs_dic_count_uc_nts(seqs_dic) |
| 237 assert c_uc_nt, "no uppercase nucleotides in input .fa sequences. "\ | 259 assert c_uc_nt, ( |
| 238 "Please change sequences to uppercase (keep in mind "\ | 260 "no uppercase nucleotides in input .fa sequences. " |
| 239 "GraphProt only scores uppercase regions (according "\ | 261 "Please change sequences to uppercase (keep in mind " |
| 240 "to its viewpoint concept)" | 262 "GraphProt only scores uppercase regions (according " |
| 263 "to its viewpoint concept)" | |
| 264 ) | |
| 241 if not args.ws_pred: | 265 if not args.ws_pred: |
| 242 # Check for lowercase sequences. | 266 # Check for lowercase sequences. |
| 243 c_lc_nt = gplib.seqs_dic_count_lc_nts(seqs_dic) | 267 c_lc_nt = gplib.seqs_dic_count_lc_nts(seqs_dic) |
| 244 assert not c_lc_nt, "lowercase nucleotides in input .fa not allowed"\ | 268 assert not c_lc_nt, ( |
| 245 " in profile predictions, since GraphProt only scores "\ | 269 "lowercase nucleotides in input .fa not allowed" |
| 270 " in profile predictions, since GraphProt only scores " | |
| 246 "uppercase regions (according to its viewpoint concept)" | 271 "uppercase regions (according to its viewpoint concept)" |
| 272 ) | |
| 247 # Check .bed. | 273 # Check .bed. |
| 248 if args.genomic_sites_bed: | 274 if args.genomic_sites_bed: |
| 249 # An array of checks, marvelous. | 275 # An array of checks, marvelous. |
| 250 assert \ | 276 assert os.path.exists( |
| 251 os.path.exists(args.genomic_sites_bed), \ | 277 args.genomic_sites_bed |
| 252 "genomic .bed file \"%s\" not found" % (args.genomic_sites_bed) | 278 ), 'genomic .bed file "%s" not found' % (args.genomic_sites_bed) |
| 253 # Check .bed for content. | 279 # Check .bed for content. |
| 254 assert gplib.count_file_rows(args.genomic_sites_bed), \ | 280 assert gplib.count_file_rows( |
| 255 "genomic .bed file \"%s\" is empty" % (args.genomic_sites_bed) | 281 args.genomic_sites_bed |
| 282 ), 'genomic .bed file "%s" is empty' % (args.genomic_sites_bed) | |
| 256 # Check .bed for 6-column format. | 283 # Check .bed for 6-column format. |
| 257 assert gplib.bed_check_six_col_format(args.genomic_sites_bed), \ | 284 assert gplib.bed_check_six_col_format( |
| 258 "genomic .bed file \"%s\" appears to not "\ | 285 args.genomic_sites_bed |
| 259 "be in 6-column .bed format" % (args.genomic_sites_bed) | 286 ), 'genomic .bed file "%s" appears to not ' "be in 6-column .bed format" % ( |
| 287 args.genomic_sites_bed | |
| 288 ) | |
| 260 # Check for unique column 4 IDs. | 289 # Check for unique column 4 IDs. |
| 261 assert gplib.bed_check_unique_ids(args.genomic_sites_bed), \ | 290 assert gplib.bed_check_unique_ids( |
| 262 "genomic .bed file \"%s\" column 4 IDs not unique" % \ | 291 args.genomic_sites_bed |
| 263 (args.genomic_sites_bed) | 292 ), 'genomic .bed file "%s" column 4 IDs not unique' % (args.genomic_sites_bed) |
| 264 # Read in .bed regions, compare to FASTA sequences. | 293 # Read in .bed regions, compare to FASTA sequences. |
| 265 seq_len_dic = gplib.get_seq_lengths_from_seqs_dic(seqs_dic) | 294 seq_len_dic = gplib.get_seq_lengths_from_seqs_dic(seqs_dic) |
| 266 reg_len_dic = gplib.bed_get_region_lengths(args.genomic_sites_bed) | 295 reg_len_dic = gplib.bed_get_region_lengths(args.genomic_sites_bed) |
| 267 for seq_id in seq_len_dic: | 296 for seq_id in seq_len_dic: |
| 268 seq_l = seq_len_dic[seq_id] | 297 seq_l = seq_len_dic[seq_id] |
| 269 assert seq_id in reg_len_dic, \ | 298 assert ( |
| 270 "sequence ID \"%s\" missing in input .bed \"%s\"" % \ | 299 seq_id in reg_len_dic |
| 271 (seq_id, args.genomic_sites_bed) | 300 ), 'sequence ID "%s" missing in input .bed "%s"' % ( |
| 301 seq_id, | |
| 302 args.genomic_sites_bed, | |
| 303 ) | |
| 272 reg_l = reg_len_dic[seq_id] | 304 reg_l = reg_len_dic[seq_id] |
| 273 assert seq_l == reg_l, \ | 305 assert ( |
| 274 "sequence length differs from .bed region length (%i != %i)" \ | 306 seq_l == reg_l |
| 275 % (seq_l, reg_l) | 307 ), "sequence length differs from .bed region length (%i != %i)" % ( |
| 308 seq_l, | |
| 309 reg_l, | |
| 310 ) | |
| 276 # Read in model parameters. | 311 # Read in model parameters. |
| 277 param_dic = gplib.graphprot_get_param_dic(args.in_params) | 312 param_dic = gplib.graphprot_get_param_dic(args.in_params) |
| 278 # Create GraphProt parameter string. | 313 # Create GraphProt parameter string. |
| 279 param_string = gplib.graphprot_get_param_string(args.in_params) | 314 param_string = gplib.graphprot_get_param_string(args.in_params) |
| 280 | 315 |
| 282 Run predictions. | 317 Run predictions. |
| 283 | 318 |
| 284 """ | 319 """ |
| 285 if args.ws_pred: | 320 if args.ws_pred: |
| 286 # Do whole site prediction. | 321 # Do whole site prediction. |
| 287 print("Starting whole site predictions on input .fa file" | 322 print( |
| 288 " (-action predict) ... ") | 323 "Starting whole site predictions on input .fa file" |
| 289 check_cmd = "GraphProt.pl -action predict -prefix " \ | 324 " (-action predict) ... " |
| 290 + args.data_id + " -fasta " + args.in_fa + " " \ | 325 ) |
| 291 + param_string + " -model " + args.in_model | 326 check_cmd = ( |
| 327 "GraphProt.pl -action predict -prefix " | |
| 328 + args.data_id | |
| 329 + " -fasta " | |
| 330 + args.in_fa | |
| 331 + " " | |
| 332 + param_string | |
| 333 + " -model " | |
| 334 + args.in_model | |
| 335 ) | |
| 292 output = subprocess.getoutput(check_cmd) | 336 output = subprocess.getoutput(check_cmd) |
| 293 assert output, "the following call of GraphProt.pl "\ | 337 assert ( |
| 294 "produced no output:\n%s" % (check_cmd) | 338 output |
| 339 ), "the following call of GraphProt.pl " "produced no output:\n%s" % (check_cmd) | |
| 295 if args.gp_output: | 340 if args.gp_output: |
| 296 print(output) | 341 print(output) |
| 297 ws_predictions_file = args.data_id + ".predictions" | 342 ws_predictions_file = args.data_id + ".predictions" |
| 298 assert os.path.exists(ws_predictions_file), \ | 343 assert os.path.exists( |
| 299 "Whole site prediction output .predictions file \"%s\" not found" \ | 344 ws_predictions_file |
| 300 % (ws_predictions_file) | 345 ), 'Whole site prediction output .predictions file "%s" not found' % ( |
| 346 ws_predictions_file | |
| 347 ) | |
| 301 if args.conf_out: | 348 if args.conf_out: |
| 302 # Filter by pos_train_ws_pred_median median. | 349 # Filter by pos_train_ws_pred_median median. |
| 303 assert "pos_train_ws_pred_median" in param_dic, \ | 350 assert "pos_train_ws_pred_median" in param_dic, ( |
| 304 "whole site top scores median information "\ | 351 "whole site top scores median information " "missing in .params file" |
| 305 "missing in .params file" | 352 ) |
| 306 pos_tr_ws_pred_med = \ | 353 pos_tr_ws_pred_med = float(param_dic["pos_train_ws_pred_median"]) |
| 307 float(param_dic["pos_train_ws_pred_median"]) | |
| 308 # Filtered file. | 354 # Filtered file. |
| 309 filt_ws_predictions_file = args.data_id + ".p50.predictions" | 355 filt_ws_predictions_file = args.data_id + ".p50.predictions" |
| 310 print("Extracting p50 sites from whole site predictions" | 356 print( |
| 311 " (score threshold = %f) ... " % (pos_tr_ws_pred_med)) | 357 "Extracting p50 sites from whole site predictions" |
| 312 gplib.graphprot_filter_predictions_file(ws_predictions_file, | 358 " (score threshold = %f) ... " % (pos_tr_ws_pred_med) |
| 313 filt_ws_predictions_file, | 359 ) |
| 314 sc_thr=pos_tr_ws_pred_med) | 360 gplib.graphprot_filter_predictions_file( |
| 361 ws_predictions_file, filt_ws_predictions_file, sc_thr=pos_tr_ws_pred_med | |
| 362 ) | |
| 315 else: | 363 else: |
| 316 # Do profile prediction. | 364 # Do profile prediction. |
| 317 print("Starting profile predictions on on input .fa file " | 365 print( |
| 318 "(-action predict_profile) ... ") | 366 "Starting profile predictions on on input .fa file " |
| 319 check_cmd = "GraphProt.pl -action predict_profile -prefix " \ | 367 "(-action predict_profile) ... " |
| 320 + args.data_id + " -fasta " + args.in_fa + " " + \ | 368 ) |
| 321 param_string + " -model " + args.in_model | 369 check_cmd = ( |
| 370 "GraphProt.pl -action predict_profile -prefix " | |
| 371 + args.data_id | |
| 372 + " -fasta " | |
| 373 + args.in_fa | |
| 374 + " " | |
| 375 + param_string | |
| 376 + " -model " | |
| 377 + args.in_model | |
| 378 ) | |
| 322 output = subprocess.getoutput(check_cmd) | 379 output = subprocess.getoutput(check_cmd) |
| 323 assert output, "the following call of GraphProt.pl produced "\ | 380 assert ( |
| 324 "no output:\n%s" % (check_cmd) | 381 output |
| 382 ), "the following call of GraphProt.pl produced " "no output:\n%s" % (check_cmd) | |
| 325 if args.gp_output: | 383 if args.gp_output: |
| 326 print(output) | 384 print(output) |
| 327 profile_predictions_file = args.data_id + ".profile" | 385 profile_predictions_file = args.data_id + ".profile" |
| 328 assert os.path.exists(profile_predictions_file), \ | 386 assert os.path.exists( |
| 329 "Profile prediction output .profile file \"%s\" not found" % \ | 387 profile_predictions_file |
| 330 (profile_predictions_file) | 388 ), 'Profile prediction output .profile file "%s" not found' % ( |
| 389 profile_predictions_file | |
| 390 ) | |
| 331 | 391 |
| 332 # Profile prediction output files. | 392 # Profile prediction output files. |
| 333 avg_prof_file = args.data_id + ".avg_profile" | 393 avg_prof_file = args.data_id + ".avg_profile" |
| 334 avg_prof_peaks_file = args.data_id + ".avg_profile.peaks.bed" | 394 avg_prof_peaks_file = args.data_id + ".avg_profile.peaks.bed" |
| 335 avg_prof_gen_peaks_file = args.data_id + \ | 395 avg_prof_gen_peaks_file = args.data_id + ".avg_profile.genomic_peaks.bed" |
| 336 ".avg_profile.genomic_peaks.bed" | 396 avg_prof_peaks_p50_file = args.data_id + ".avg_profile.p50.peaks.bed" |
| 337 avg_prof_peaks_p50_file = args.data_id + \ | 397 avg_prof_gen_peaks_p50_file = ( |
| 338 ".avg_profile.p50.peaks.bed" | 398 args.data_id + ".avg_profile.p50.genomic_peaks.bed" |
| 339 avg_prof_gen_peaks_p50_file = args.data_id + \ | 399 ) |
| 340 ".avg_profile.p50.genomic_peaks.bed" | |
| 341 | 400 |
| 342 # Get sequence IDs in order from input .fa file. | 401 # Get sequence IDs in order from input .fa file. |
| 343 seq_ids_list = gplib.fasta_read_in_ids(args.in_fa) | 402 seq_ids_list = gplib.fasta_read_in_ids(args.in_fa) |
| 344 # Calculate average profiles. | 403 # Calculate average profiles. |
| 345 print("Getting average profile from profile " | 404 print( |
| 346 "(extlr for smoothing: %i) ... " % (args.ap_extlr)) | 405 "Getting average profile from profile " |
| 347 gplib.graphprot_profile_calc_avg_profile(profile_predictions_file, | 406 "(extlr for smoothing: %i) ... " % (args.ap_extlr) |
| 348 avg_prof_file, | 407 ) |
| 349 ap_extlr=args.ap_extlr, | 408 gplib.graphprot_profile_calc_avg_profile( |
| 350 seq_ids_list=seq_ids_list, | 409 profile_predictions_file, |
| 351 method=2) | 410 avg_prof_file, |
| 411 ap_extlr=args.ap_extlr, | |
| 412 seq_ids_list=seq_ids_list, | |
| 413 method=2, | |
| 414 ) | |
| 352 # Extract peak regions on sequences with threshold score 0. | 415 # Extract peak regions on sequences with threshold score 0. |
| 353 print("Extracting peak regions from average profile " | 416 print( |
| 354 "(score threshold = 0) ... ") | 417 "Extracting peak regions from average profile " "(score threshold = 0) ... " |
| 418 ) | |
| 355 killpep8 = args.max_merge_dist | 419 killpep8 = args.max_merge_dist |
| 356 gplib.graphprot_profile_extract_peak_regions(avg_prof_file, | 420 gplib.graphprot_profile_extract_peak_regions( |
| 357 avg_prof_peaks_file, | 421 avg_prof_file, |
| 358 max_merge_dist=killpep8, | 422 avg_prof_peaks_file, |
| 359 sc_thr=args.score_thr) | 423 max_merge_dist=killpep8, |
| 424 sc_thr=args.score_thr, | |
| 425 ) | |
| 360 # Convert peaks to genomic coordinates. | 426 # Convert peaks to genomic coordinates. |
| 361 if args.genomic_sites_bed: | 427 if args.genomic_sites_bed: |
| 362 print("Converting peak regions to genomic coordinates ... ") | 428 print("Converting peak regions to genomic coordinates ... ") |
| 363 killit = args.genomic_sites_bed | 429 killit = args.genomic_sites_bed |
| 364 gplib.bed_peaks_to_genomic_peaks(avg_prof_peaks_file, | 430 gplib.bed_peaks_to_genomic_peaks( |
| 365 avg_prof_gen_peaks_file, | 431 avg_prof_peaks_file, |
| 366 print_rows=False, | 432 avg_prof_gen_peaks_file, |
| 367 genomic_sites_bed=killit) | 433 print_rows=False, |
| 434 genomic_sites_bed=killit, | |
| 435 ) | |
| 368 # Extract peak regions with threshold score p50. | 436 # Extract peak regions with threshold score p50. |
| 369 if args.conf_out: | 437 if args.conf_out: |
| 370 sc_id = "pos_train_avg_profile_median_%i" % (args.ap_extlr) | 438 sc_id = "pos_train_avg_profile_median_%i" % (args.ap_extlr) |
| 371 # Filter by pos_tr_ws_pred_med median. | 439 # Filter by pos_tr_ws_pred_med median. |
| 372 assert sc_id in param_dic, "average profile extlr %i median "\ | 440 assert sc_id in param_dic, ( |
| 441 "average profile extlr %i median " | |
| 373 "information missing in .params file" % (args.ap_extlr) | 442 "information missing in .params file" % (args.ap_extlr) |
| 443 ) | |
| 374 p50_sc_thr = float(param_dic[sc_id]) | 444 p50_sc_thr = float(param_dic[sc_id]) |
| 375 print("Extracting p50 peak regions from average profile " | 445 print( |
| 376 "(score threshold = %f) ... " % (p50_sc_thr)) | 446 "Extracting p50 peak regions from average profile " |
| 447 "(score threshold = %f) ... " % (p50_sc_thr) | |
| 448 ) | |
| 377 despair = avg_prof_peaks_p50_file | 449 despair = avg_prof_peaks_p50_file |
| 378 pain = args.max_merge_dist | 450 pain = args.max_merge_dist |
| 379 gplib.graphprot_profile_extract_peak_regions(avg_prof_file, | 451 gplib.graphprot_profile_extract_peak_regions( |
| 380 despair, | 452 avg_prof_file, despair, max_merge_dist=pain, sc_thr=p50_sc_thr |
| 381 max_merge_dist=pain, | 453 ) |
| 382 sc_thr=p50_sc_thr) | |
| 383 # Convert peaks to genomic coordinates. | 454 # Convert peaks to genomic coordinates. |
| 384 if args.genomic_sites_bed: | 455 if args.genomic_sites_bed: |
| 385 print("Converting p50 peak regions to " | 456 print("Converting p50 peak regions to " "genomic coordinates ... ") |
| 386 "genomic coordinates ... ") | |
| 387 madness = args.genomic_sites_bed | 457 madness = args.genomic_sites_bed |
| 388 gplib.bed_peaks_to_genomic_peaks(avg_prof_peaks_p50_file, | 458 gplib.bed_peaks_to_genomic_peaks( |
| 389 avg_prof_gen_peaks_p50_file, | 459 avg_prof_peaks_p50_file, |
| 390 genomic_sites_bed=madness) | 460 avg_prof_gen_peaks_p50_file, |
| 461 genomic_sites_bed=madness, | |
| 462 ) | |
| 391 # Done. | 463 # Done. |
| 392 print("Script: I'm done.") | 464 print("Script: I'm done.") |
| 393 print("Author: ... ") | 465 print("Author: ... ") |