Merge pull request #46 from hiruna72/dev

Dev
hiruna72 · Oct 6, 2023 · 42643f5 · 42643f5
2 parents 1110964 + f9674c5
commit 42643f5
Show file tree

Hide file tree

Showing 6 changed files with 54 additions and 9 deletions.
diff --git a/README.md b/README.md
@@ -31,8 +31,9 @@ Google Chrome is the recommended web browser to visualise these plots.
    1. [Base shift](#base-shift)
    2. [Signal scaling](#signal-scaling)
 11. [Plot conventions](#plot-conventions)
-12. [Notes](#notes)
-13. [Examples](#examples)
+12. [Calculate alignment statistics](#calculate-alignment-statistics)
+13. [Notes](#notes)
+14. [Examples](#examples)
 
 
 ## Quickstart
@@ -545,6 +546,10 @@ The commands `plot` and `plot_pileup` can take the argument `--sig_scale`. By pr
 * **D** indicates the true sequencing direction of the signals.
 
 
+## Calculate alignment statistics
+Calculate basic statistics of read/reference - signal alignments.
+Check [here](docs/commands.md/#metric) for the command.
+Check [here](docs/different_alignments.md) for an example.
 
 ## Notes
 

diff --git a/docs/commands.md b/docs/commands.md
@@ -114,6 +114,8 @@ Plot read/reference - signal alignments.
 	The implementation of each method can be found at `src/plot_utils.py/scale_signal()`
 * `--no_pa`
 	(optional) Do not convert the signal to pA levels. By default, the raw signal is converted to pA levels [default value: false].
+* `--loose_bound`
+	(optional) Also plot alignments not completely within the specified region but at least part is [default value: false].
 * `--point_size INT`
 	(optional) Radius of the signal point drawn in the plot [default value: 0.5].
 * `--base_width INT`
@@ -301,6 +303,8 @@ Instead of generating figures `metric` will generate statistics after parsing th
 	The implementation of each method can be found at `src/plot_utils.py/scale_signal()`
 * `--no_pa`
 	(optional) Do not convert the signal to pA levels. By default, the raw signal is converted to pA levels [default value: false].
+* `--loose_bound`
+	(optional) Also plot alignments not completely within the specified region but at least part is [default value: false].
 * `--base_shift INT`
 	(optional) The number of bases to shift to align fist signal move [default value: 0]. More information on this can be found at [here](pore_model.md)
 *  `--profile STR`:<br/>

diff --git a/src/_version.py b/src/_version.py
@@ -1,4 +1,4 @@
 """version.py
 Version of the software
 """
-__version__ = "0.3.0-beta"
+__version__ = "0.4.0"
diff --git a/src/metric.py b/src/metric.py
@@ -489,6 +489,15 @@ def run(args):
             else:
                 base_limit = BASE_LIMIT
             sam_record_reference_end = reference_start + ref_seq_len #1based closed
+            if not args.loose_bound:
+                if data_is_rna == 1:
+                    if args_ref_start < reference_start + 1 - kmer_correction:
+                        continue
+                else:
+                    if args_ref_start < reference_start + 1:
+                        continue
+                if args_ref_end > sam_record_reference_end:
+                    continue
 
             ref_name = sam_record.reference_name
             ref_start = reference_start + 1
@@ -685,7 +694,16 @@ def run(args):
             else:
                 base_limit = BASE_LIMIT
             paf_record_reference_end = reference_start + ref_seq_len #1based closed
-
+            if not args.loose_bound:
+                if data_is_rna == 1:
+                    if args_ref_start < reference_start + 1 - kmer_correction:
+                        continue
+                else:
+                    if args_ref_start < reference_start + 1:
+                        continue
+                if args_ref_end > paf_record_reference_end:
+                    continue
+
             ref_name = paf_record[SEQUENCE_ID]
             ref_start = reference_start + 1
             ref_end = ref_start + base_limit - 1 #ref_end is 1based closed
@@ -849,6 +867,7 @@ def argparser():
     parser.add_argument('--sig_scale', required=False, type=str, default="", help="plot the scaled signal. Supported scalings: [medmad, znorm, scaledpA]")
     # parser.add_argument('--reverse_signal', required=False, action='store_true', help="plot RNA reference/read from 5`-3` and reverse the signal")
     parser.add_argument('--no_pa', required=False, action='store_false', help="skip converting the signal to pA values")
+    parser.add_argument('--loose_bound', required=False, action='store_true', help="also plot alignments not completely within the specified region")
     parser.add_argument('--base_shift', required=False, type=int, default=PLOT_BASE_SHIFT, help="the number of bases to shift to align fist signal move")
     parser.add_argument('--profile', required=False, default="", type=str, help="determine base_shift using preset values")
     parser.add_argument('--list_profile', action='store_true', help="list the available profiles")

diff --git a/src/plot.py b/src/plot.py
@@ -776,6 +776,15 @@ def run(args):
             else:
                 base_limit = BASE_LIMIT
             sam_record_reference_end = reference_start + ref_seq_len #1based closed
+            if not args.loose_bound:
+                if data_is_rna == 1:
+                    if args_ref_start < reference_start + 1 - kmer_correction:
+                        continue
+                else:
+                    if args_ref_start < reference_start + 1:
+                        continue
+                if args_ref_end > sam_record_reference_end:
+                    continue
 
             ref_name = sam_record.reference_name
             ref_start = reference_start + 1
@@ -984,6 +993,15 @@ def run(args):
             else:
                 base_limit = BASE_LIMIT
             paf_record_reference_end = reference_start + ref_seq_len #1based closed
+            if not args.loose_bound:
+                if data_is_rna == 1:
+                    if args_ref_start < reference_start + 1 - kmer_correction:
+                        continue
+                else:
+                    if args_ref_start < reference_start + 1:
+                        continue
+                if args_ref_end > paf_record_reference_end:
+                    continue
 
             ref_name = paf_record[SEQUENCE_ID]
             ref_start = reference_start + 1
@@ -1160,6 +1178,7 @@ def argparser():
     parser.add_argument('--sig_scale', required=False, type=str, default="", help="plot the scaled signal. Supported scalings: [medmad, znorm, scaledpA]")
     # parser.add_argument('--reverse_signal', required=False, action='store_true', help="plot RNA reference/read from 5`-3` and reverse the signal")
     parser.add_argument('--no_pa', required=False, action='store_false', help="skip converting the signal to pA values")
+    parser.add_argument('--loose_bound', required=False, action='store_true', help="also plot alignments not completely within the specified region")
     parser.add_argument('--point_size', required=False, type=int, default=0.5, help="signal point radius [0.5]")
     parser.add_argument('--base_width', required=False, type=int, default=FIXED_BASE_WIDTH, help="base width when plotting with fixed base width")
     parser.add_argument('--base_shift', required=False, type=int, default=PLOT_BASE_SHIFT, help="the number of bases to shift to align fist signal move")

diff --git a/src/plot_pileup.py b/src/plot_pileup.py
@@ -352,8 +352,6 @@ def run(args):
     else:
         if not os.path.exists(args.output_dir):
             os.mkdir(args.output_dir)
-    # if not args.loose_bound:
-    #     print("Warning: the alignments that are not completely within the specified region will be skipped. Please specify --loose_bound to keep")
 
     # open signal file
     s5 = pyslow5.Open(args.slow5, 'r')
@@ -456,7 +454,7 @@ def run(args):
             else:
                 base_limit = BASE_LIMIT
             sam_record_reference_end = reference_start + ref_seq_len #1based closed
-            if not args.loose_bound:
+            if True: # if not args.loose_bound:
                 if data_is_rna == 1:
                     if args_ref_start < reference_start + 1 - kmer_correction:
                         continue
@@ -690,7 +688,7 @@ def run(args):
             else:
                 base_limit = BASE_LIMIT
             paf_record_reference_end = reference_start + ref_seq_len #1based closed
-            if not args.loose_bound:
+            if True: # if not args.loose_bound:
                 if data_is_rna == 1:
                     if args_ref_start < reference_start + 1 - kmer_correction:
                         continue
@@ -948,7 +946,7 @@ def argparser():
     parser.add_argument('--overlap_bottom', required=False, action='store_true', help="plot the overlap at the bottom")
     parser.add_argument('--no_overlap', required=False, action='store_true', help="skip plotting the overlap")
     parser.add_argument('--overlap_only', required=False, action='store_true', help="plot only the overlap")
-    parser.add_argument('--loose_bound', required=False, action='store_true', help="also plot alignments not completely within the specified region")
+    # parser.add_argument('--loose_bound', required=False, action='store_true', help="also plot alignments not completely within the specified region")
     parser.add_argument('--point_size', required=False, type=int, default=0.5, help="signal point radius [0.5]")
     parser.add_argument('--base_width', required=False, type=int, default=FIXED_BASE_WIDTH, help="base width when plotting with fixed base width")
     parser.add_argument('--base_shift', required=False, type=int, default=PLOT_BASE_SHIFT, help="the number of bases to shift to align fist signal move")