diff gplib.py @ 5:ddcf35a868b8 draft

planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/tools/rna_tools/graphprot commit ad60258f5759eaa205fec4af6143c728ea131419
author bgruening
date Wed, 05 Jun 2024 16:40:51 +0000
parents ace92c9a4653
children
line wrap: on
line diff
--- a/gplib.py	Thu Jan 28 15:06:14 2021 +0000
+++ b/gplib.py	Wed Jun 05 16:40:51 2024 +0000
@@ -1,4 +1,3 @@
-
 import gzip
 import random
 import re
@@ -6,7 +5,6 @@
 import subprocess
 from distutils.spawn import find_executable
 
-
 """
 
 Run doctests:
@@ -17,7 +15,8 @@
 """
 
 
-###############################################################################
+#######################################################################
+
 
 def graphprot_predictions_get_median(predictions_file):
     """
@@ -41,11 +40,12 @@
     return statistics.median(sc_list)
 
 
-###############################################################################
+#######################################################################
+
 
-def graphprot_profile_get_tsm(profile_file,
-                              profile_type="profile",
-                              avg_profile_extlr=5):
+def graphprot_profile_get_tsm(
+    profile_file, profile_type="profile", avg_profile_extlr=5
+):
 
     """
     Given a GraphProt .profile file, extract for each site (identified by
@@ -88,23 +88,21 @@
             max_list.append(max_sc)
         elif profile_type == "avg_profile":
             # Convert profile score list to average profile scores list.
-            aps_list = \
-                list_moving_window_average_values(lists_dic[seq_id],
-                                                  win_extlr=avg_profile_extlr)
+            aps_list = list_moving_window_average_values(
+                lists_dic[seq_id], win_extlr=avg_profile_extlr
+            )
             max_sc = max(aps_list)
             max_list.append(max_sc)
         else:
-            assert 0, "invalid profile_type argument given: \"%s\"" \
-                % (profile_type)
+            assert 0, 'invalid profile_type argument given: "%s"' % (profile_type)
     # Return the median.
     return statistics.median(max_list)
 
 
-###############################################################################
+#######################################################################
 
-def list_moving_window_average_values(in_list,
-                                      win_extlr=5,
-                                      method=1):
+
+def list_moving_window_average_values(in_list, win_extlr=5, method=1):
     """
     Take a list of numeric values, and calculate for each position a new value,
     by taking the mean value of the window of positions -win_extlr and
@@ -152,7 +150,8 @@
     return new_list
 
 
-###############################################################################
+#######################################################################
+
 
 def echo_add_to_file(echo_string, out_file):
     """
@@ -167,14 +166,16 @@
     assert not error, "echo is complaining:\n%s\n%s" % (check_cmd, output)
 
 
-###############################################################################
+#######################################################################
+
 
 def is_tool(name):
     """Check whether tool "name" is in PATH."""
     return find_executable(name) is not None
 
 
-###############################################################################
+#######################################################################
+
 
 def count_fasta_headers(fasta_file):
     """
@@ -194,7 +195,8 @@
     return row_count
 
 
-###############################################################################
+#######################################################################
+
 
 def make_file_copy(in_file, out_file):
     """
@@ -202,21 +204,26 @@
 
     """
     check_cmd = "cat " + in_file + " > " + out_file
-    assert in_file != out_file, \
-        "cat does not like to cat file into same file (%s)" % (check_cmd)
+    assert in_file != out_file, "cat does not like to cat file into same file (%s)" % (
+        check_cmd
+    )
     output = subprocess.getoutput(check_cmd)
     error = False
     if output:
         error = True
-    assert not error, \
-        "cat did not like your input (in_file: %s, out_file: %s):\n%s" \
-        % (in_file, out_file, output)
+    assert not error, "cat did not like your input (in_file: %s, out_file: %s):\n%s" % (
+        in_file,
+        out_file,
+        output,
+    )
 
 
-###############################################################################
+#######################################################################
+
 
-def split_fasta_into_test_train_files(in_fasta, test_out_fa, train_out_fa,
-                                      test_size=500):
+def split_fasta_into_test_train_files(
+    in_fasta, test_out_fa, train_out_fa, test_size=500
+):
     """
     Split in_fasta .fa file into two files (e.g. test, train).
 
@@ -230,7 +237,7 @@
     TRAINOUT = open(train_out_fa, "w")
     for seq_id in rand_ids_list:
         seq = seqs_dic[seq_id]
-        if (c_out >= test_size):
+        if c_out >= test_size:
             TRAINOUT.write(">%s\n%s\n" % (seq_id, seq))
         else:
             TESTOUT.write(">%s\n%s\n" % (seq_id, seq))
@@ -239,7 +246,8 @@
     TRAINOUT.close()
 
 
-###############################################################################
+#######################################################################
+
 
 def check_seqs_dic_format(seqs_dic):
     """
@@ -267,16 +275,19 @@
     return bad_seq_ids
 
 
-###############################################################################
+#######################################################################
+
 
-def read_fasta_into_dic(fasta_file,
-                        seqs_dic=False,
-                        ids_dic=False,
-                        read_dna=False,
-                        short_ensembl=False,
-                        reject_lc=False,
-                        convert_to_uc=False,
-                        skip_n_seqs=True):
+def read_fasta_into_dic(
+    fasta_file,
+    seqs_dic=False,
+    ids_dic=False,
+    read_dna=False,
+    short_ensembl=False,
+    reject_lc=False,
+    convert_to_uc=False,
+    skip_n_seqs=True,
+):
     """
     Read in FASTA sequences, convert to RNA, store in dictionary
     and return dictionary.
@@ -302,7 +313,7 @@
 
     # Go through FASTA file, extract sequences.
     if re.search(r".+\.gz$", fasta_file):
-        f = gzip.open(fasta_file, 'rt')
+        f = gzip.open(fasta_file, "rt")
     else:
         f = open(fasta_file, "r")
     for line in f:
@@ -315,9 +326,10 @@
                 if re.search(r".+\..+", seq_id):
                     m = re.search(r"(.+?)\..+", seq_id)
                     seq_id = m.group(1)
-            assert seq_id not in seqs_dic, \
-                "non-unique FASTA header \"%s\" in \"%s\"" \
-                % (seq_id, fasta_file)
+            assert seq_id not in seqs_dic, 'non-unique FASTA header "%s" in "%s"' % (
+                seq_id,
+                fasta_file,
+            )
             if ids_dic:
                 if seq_id in ids_dic:
                     seqs_dic[seq_id] = ""
@@ -328,31 +340,31 @@
                 m = re.search("([ACGTUN]+)", line, re.I)
                 seq = m.group(1)
                 if reject_lc:
-                    assert \
-                        not re.search("[a-z]", seq), \
-                        "lc char detected in seq \"%i\" (reject_lc=True)" \
-                        % (seq_id)
+                    assert not re.search(
+                        "[a-z]", seq
+                    ), 'lc char detected in seq "%i" (reject_lc=True)' % (seq_id)
                 if convert_to_uc:
                     seq = seq.upper()
                 # If sequences with N nucleotides should be skipped.
                 if skip_n_seqs:
                     if "n" in m.group(1) or "N" in m.group(1):
-                        print("WARNING: \"%s\" contains N. Discarding "
-                              "sequence ... " % (seq_id))
+                        print(
+                            'WARNING: "%s" contains N. Discarding '
+                            "sequence ... " % (seq_id)
+                        )
                         del seqs_dic[seq_id]
                         continue
                 # Convert to RNA, concatenate sequence.
                 if read_dna:
-                    seqs_dic[seq_id] += \
-                        m.group(1).replace("U", "T").replace("u", "t")
+                    seqs_dic[seq_id] += m.group(1).replace("U", "T").replace("u", "t")
                 else:
-                    seqs_dic[seq_id] += \
-                        m.group(1).replace("T", "U").replace("t", "u")
+                    seqs_dic[seq_id] += m.group(1).replace("T", "U").replace("t", "u")
     f.close()
     return seqs_dic
 
 
-###############################################################################
+#######################################################################
+
 
 def random_order_dic_keys_into_list(in_dic):
     """
@@ -366,7 +378,8 @@
     return id_list
 
 
-###############################################################################
+#######################################################################
+
 
 def graphprot_get_param_string(params_file):
     """
@@ -394,11 +407,12 @@
                 else:
                     param_string += "-%s %s " % (par, setting)
             else:
-                assert 0, "pattern matching failed for string \"%s\"" % (param)
+                assert 0, 'pattern matching failed for string "%s"' % (param)
     return param_string
 
 
-###############################################################################
+#######################################################################
+
 
 def seqs_dic_count_uc_nts(seqs_dic):
     """
@@ -416,11 +430,12 @@
     assert seqs_dic, "Given sequence dictionary empty"
     c_uc = 0
     for seq_id in seqs_dic:
-        c_uc += len(re.findall(r'[A-Z]', seqs_dic[seq_id]))
+        c_uc += len(re.findall(r"[A-Z]", seqs_dic[seq_id]))
     return c_uc
 
 
-###############################################################################
+#######################################################################
+
 
 def seqs_dic_count_lc_nts(seqs_dic):
     """
@@ -438,11 +453,12 @@
     assert seqs_dic, "Given sequence dictionary empty"
     c_uc = 0
     for seq_id in seqs_dic:
-        c_uc += len(re.findall(r'[a-z]', seqs_dic[seq_id]))
+        c_uc += len(re.findall(r"[a-z]", seqs_dic[seq_id]))
     return c_uc
 
 
-###############################################################################
+#######################################################################
+
 
 def count_file_rows(in_file):
     """
@@ -462,7 +478,8 @@
     return row_count
 
 
-###############################################################################
+#######################################################################
+
 
 def bed_check_six_col_format(bed_file):
     """
@@ -488,7 +505,8 @@
     return six_col_format
 
 
-###############################################################################
+#######################################################################
+
 
 def bed_check_unique_ids(bed_file):
     """
@@ -512,7 +530,8 @@
         return True
 
 
-###############################################################################
+#######################################################################
+
 
 def get_seq_lengths_from_seqs_dic(seqs_dic):
     """
@@ -527,7 +546,8 @@
     return seq_len_dic
 
 
-###############################################################################
+#######################################################################
+
 
 def bed_get_region_lengths(bed_file):
     """
@@ -548,19 +568,19 @@
             site_e = int(cols[2])
             site_id = cols[3]
             site_l = site_e - site_s
-            assert site_id \
-                not in id2len_dic, \
-                "column 4 IDs not unique in given .bed file \"%s\"" \
-                % (bed_file)
+            assert (
+                site_id not in id2len_dic
+            ), 'column 4 IDs not unique in given .bed file "%s"' % (bed_file)
             id2len_dic[site_id] = site_l
     f.closed
-    assert id2len_dic, \
-        "No IDs read into dic (input file \"%s\" empty or malformatted?)" \
-        % (bed_file)
+    assert (
+        id2len_dic
+    ), 'No IDs read into dic (input file "%s" empty or malformatted?)' % (bed_file)
     return id2len_dic
 
 
-###############################################################################
+#######################################################################
+
 
 def graphprot_get_param_dic(params_file):
     """
@@ -599,10 +619,10 @@
     return param_dic
 
 
-###############################################################################
+#######################################################################
 
-def graphprot_filter_predictions_file(in_file, out_file,
-                                      sc_thr=0):
+
+def graphprot_filter_predictions_file(in_file, out_file, sc_thr=0):
     """
     Filter GraphProt .predictions file by given score thr_sc.
     """
@@ -619,7 +639,8 @@
     OUTPRED.close()
 
 
-###############################################################################
+#######################################################################
+
 
 def fasta_read_in_ids(fasta_file):
     """
@@ -642,12 +663,12 @@
     return ids_list
 
 
-###############################################################################
+#######################################################################
+
 
-def graphprot_profile_calc_avg_profile(in_file, out_file,
-                                       ap_extlr=5,
-                                       seq_ids_list=False,
-                                       method=1):
+def graphprot_profile_calc_avg_profile(
+    in_file, out_file, ap_extlr=5, seq_ids_list=False, method=1
+):
     """
     Given a GraphProt .profile file, calculate average profiles and output
     average profile file.
@@ -702,15 +723,16 @@
         if seq_ids_list:
             c_seq_ids = len(seq_ids_list)
             c_site_ids = len(site_starts_dic)
-            assert c_seq_ids == c_site_ids, \
-                "# sequence IDs != # site IDs (%i != %i)" \
-                % (c_seq_ids, c_site_ids)
+            assert (
+                c_seq_ids == c_site_ids
+            ), "# sequence IDs != # site IDs (%i != %i)" % (c_seq_ids, c_site_ids)
         OUTPROF = open(out_file, "w")
         # For each site, calculate average profile scores list.
         for site_id in lists_dic:
             # Convert profile score list to average profile scores list.
-            aps_list = list_moving_window_average_values(lists_dic[site_id],
-                                                         win_extlr=ap_extlr)
+            aps_list = list_moving_window_average_values(
+                lists_dic[site_id], win_extlr=ap_extlr
+            )
             start_pos = site_starts_dic[site_id]
             # Get original FASTA sequence ID.
             if seq_ids_list:
@@ -741,10 +763,9 @@
                 if cur_id != old_id:
                     # Process old id scores.
                     if scores_list:
-                        aps_list = \
-                            list_moving_window_average_values(
-                                scores_list,
-                                win_extlr=ap_extlr)
+                        aps_list = list_moving_window_average_values(
+                            scores_list, win_extlr=ap_extlr
+                        )
                         start_pos = site_starts_dic[old_id]
                         seq_id = old_id
                         # Get original FASTA sequence ID.
@@ -763,8 +784,9 @@
         f.close()
         # Process last block.
         if scores_list:
-            aps_list = list_moving_window_average_values(scores_list,
-                                                         win_extlr=ap_extlr)
+            aps_list = list_moving_window_average_values(
+                scores_list, win_extlr=ap_extlr
+            )
             start_pos = site_starts_dic[old_id]
             seq_id = old_id
             # Get original FASTA sequence ID.
@@ -776,11 +798,12 @@
         OUTPROF.close()
 
 
-###############################################################################
+#######################################################################
+
 
-def graphprot_profile_extract_peak_regions(in_file, out_file,
-                                           max_merge_dist=0,
-                                           sc_thr=0):
+def graphprot_profile_extract_peak_regions(
+    in_file, out_file, max_merge_dist=0, sc_thr=0
+):
     """
     Extract peak regions from GraphProt .profile file.
     Store the peak regions (defined as regions with scores >= sc_thr)
@@ -835,21 +858,22 @@
                 # Process old id scores.
                 if scores_list:
                     # Extract peaks from region.
-                    peak_list = \
-                        list_extract_peaks(scores_list,
-                                           max_merge_dist=max_merge_dist,
-                                           coords="bed",
-                                           sc_thr=sc_thr)
+                    peak_list = list_extract_peaks(
+                        scores_list,
+                        max_merge_dist=max_merge_dist,
+                        coords="bed",
+                        sc_thr=sc_thr,
+                    )
                     start_pos = site_starts_dic[old_id]
                     # Print out peaks in .bed format.
                     for ln in peak_list:
                         peak_s = start_pos + ln[0]
                         peak_e = start_pos + ln[1]
                         site_id = "%s,%i" % (old_id, ln[2])
-                        OUTPEAKS.write("%s\t%i\t%i"
-                                       "\t%s\t%f\t+\n"
-                                       % (old_id, peak_s,
-                                          peak_e, site_id, ln[3]))
+                        OUTPEAKS.write(
+                            "%s\t%i\t%i"
+                            "\t%s\t%f\t+\n" % (old_id, peak_s, peak_e, site_id, ln[3])
+                        )
                     # Reset list.
                     scores_list = []
                 old_id = cur_id
@@ -861,27 +885,25 @@
     # Process last block.
     if scores_list:
         # Extract peaks from region.
-        peak_list = list_extract_peaks(scores_list,
-                                       max_merge_dist=max_merge_dist,
-                                       coords="bed",
-                                       sc_thr=sc_thr)
+        peak_list = list_extract_peaks(
+            scores_list, max_merge_dist=max_merge_dist, coords="bed", sc_thr=sc_thr
+        )
         start_pos = site_starts_dic[old_id]
         # Print out peaks in .bed format.
         for ln in peak_list:
             peak_s = start_pos + ln[0]
             peak_e = start_pos + ln[1]
             site_id = "%s,%i" % (old_id, ln[2])  # best score also 1-based.
-            OUTPEAKS.write("%s\t%i\t%i\t%s\t%f\t+\n"
-                           % (old_id, peak_s, peak_e, site_id, ln[3]))
+            OUTPEAKS.write(
+                "%s\t%i\t%i\t%s\t%f\t+\n" % (old_id, peak_s, peak_e, site_id, ln[3])
+            )
     OUTPEAKS.close()
 
 
-###############################################################################
+#######################################################################
 
-def list_extract_peaks(in_list,
-                       max_merge_dist=0,
-                       coords="list",
-                       sc_thr=0):
+
+def list_extract_peaks(in_list, max_merge_dist=0, coords="list", sc_thr=0):
     """
     Extract peak regions from list.
     Peak region is defined as region >= score threshold.
@@ -969,8 +991,12 @@
                     if peak_list[i][3] < peak_list[j][3]:
                         new_top_pos = peak_list[j][2]
                         new_top_sc = peak_list[j][3]
-                    new_peak = [peak_list[i][0], peak_list[j][1],
-                                new_top_pos, new_top_sc]
+                    new_peak = [
+                        peak_list[i][0],
+                        peak_list[j][1],
+                        new_top_pos,
+                        new_top_sc,
+                    ]
                 # If two peaks were merged.
                 if new_peak:
                     merged_peak_list.append(new_peak)
@@ -991,10 +1017,12 @@
     return peak_list
 
 
-###############################################################################
+#######################################################################
+
 
-def bed_peaks_to_genomic_peaks(peak_file, genomic_peak_file, genomic_sites_bed,
-                               print_rows=False):
+def bed_peaks_to_genomic_peaks(
+    peak_file, genomic_peak_file, genomic_sites_bed, print_rows=False
+):
     """
     Given a .bed file of sequence peak regions (possible coordinates from
     0 to length of s), convert peak coordinates to genomic coordinates.
@@ -1017,10 +1045,9 @@
             row = line.strip()
             cols = line.strip().split("\t")
             site_id = cols[3]
-            assert site_id \
-                not in id2row_dic, \
-                "column 4 IDs not unique in given .bed file \"%s\"" \
-                % (genomic_sites_bed)
+            assert (
+                site_id not in id2row_dic
+            ), 'column 4 IDs not unique in given .bed file "%s"' % (genomic_sites_bed)
             id2row_dic[site_id] = row
     f.close()
 
@@ -1034,13 +1061,14 @@
             site_e = int(cols[2])
             site_id2 = cols[3]
             site_sc = float(cols[4])
-            assert re.search(".+,.+", site_id2), \
-                "regular expression failed for ID \"%s\"" % (site_id2)
+            assert re.search(
+                ".+,.+", site_id2
+            ), 'regular expression failed for ID "%s"' % (site_id2)
             m = re.search(r".+,(\d+)", site_id2)
             sc_pos = int(m.group(1))  # 1-based.
-            assert site_id in id2row_dic, \
-                "site ID \"%s\" not found in genomic sites dictionary" \
-                % (site_id)
+            assert (
+                site_id in id2row_dic
+            ), 'site ID "%s" not found in genomic sites dictionary' % (site_id)
             row = id2row_dic[site_id]
             rowl = row.split("\t")
             gen_chr = rowl[0]
@@ -1054,16 +1082,23 @@
                 new_s = gen_e - site_e
                 new_e = gen_e - site_s
                 new_sc_pos = gen_e - sc_pos + 1  # keep 1-based.
-            new_row = "%s\t%i\t%i\t%s,%i\t%f\t%s" \
-                      % (gen_chr, new_s, new_e,
-                         site_id, new_sc_pos, site_sc, gen_pol)
+            new_row = "%s\t%i\t%i\t%s,%i\t%f\t%s" % (
+                gen_chr,
+                new_s,
+                new_e,
+                site_id,
+                new_sc_pos,
+                site_sc,
+                gen_pol,
+            )
             OUTPEAKS.write("%s\n" % (new_row))
             if print_rows:
                 print(new_row)
     OUTPEAKS.close()
 
 
-###############################################################################
+#######################################################################
+
 
 def diff_two_files_identical(file1, file2):
     """
@@ -1087,4 +1122,4 @@
     return same
 
 
-###############################################################################
+#######################################################################