using readpaf locally

Author: hiruna72

src/calculate_offsets.py

@@ -1,6 +1,6 @@
 import os
 import argparse
-from readpaf import parse_paf
+from src import readpaf_local
 import pysam
 import pyslow5
 from pyfaidx import Fasta
@@ -136,7 +136,7 @@ def calculate_offsets(args, sig_move_offset, output_pdf, s5):
         else:
             sequence_reads = Fastq(sequence_file)
         num_reads = 0
-        for paf_record in parse_paf(handle):
+        for paf_record in readpaf_local.parse_paf(handle):
             read_id = paf_record.query_name
             if args.read_id != "" and args.read_id != read_id:
                 continue

src/metric.py

@@ -7,7 +7,7 @@
 import pyslow5
 import argparse
 import re
-from readpaf import parse_paf
+from src import readpaf_local
 from pyfaidx import Fasta
 from pyfastx import Fastq
 import os
@@ -292,7 +292,7 @@ def run(args):
                 sequence_reads = Fasta(args.file)
             else:
                 sequence_reads = Fastq(args.file)
-            for paf_record in parse_paf(handle):
+            for paf_record in readpaf_local.parse_paf(handle):
                 metric_record = {}
                 if paf_record.query_name != paf_record.target_name:
                     raise Exception("Error: this paf file is a signal to reference mapping. Please provide the argument --sig_ref ")
@@ -880,19 +880,28 @@ def run(args):
     else:
         raise Exception("Error: You should not have ended up here. Please check your arguments")
 
+    flag_discard = 0
     if args.output_current_column:
         for key, item in column_raw_samples.items():
             # print(', '.join(map(str, item)))
             array = np.array(item)
 
             median = np.median(array)
+            if median == 0:
+                print("median is zero!")
+                flag_discard = 1
+                break
             darray = array/median
+
             print(', '.join(map(str, darray)))
 
             # darray = plot_utils.scale_signal(array, 'medmad', {})
             # print(', '.join(map(str, darray)))
 
-    print("Number of records: {}".format(num_plots))
+    if flag_discard:
+        print("discard")
+    else:
+        print("Number of records: {}".format(num_plots))
     if num_plots == 0:
         print("Squigualiser only plots reads that span across the specified region entirely. Reduce the region interval and double check with IGV : {}".format(num_plots))
 

src/plot.py

@@ -14,7 +14,7 @@
 import argparse
 import re
 import logging
-from readpaf import parse_paf
+from src import readpaf_local
 from sys import stdin
 from pyfaidx import Fasta
 from pyfastx import Fastq
@@ -669,7 +669,7 @@ def run(args):
                 args_ref_start = None
                 args_ref_end = None
 
-            for paf_record in parse_paf(handle):
+            for paf_record in readpaf_local.parse_paf(handle):
                 if args.auto:
                     kmer_correction = 0
                     draw_data["base_shift"] = 0

src/plot_pileup.py

@@ -15,7 +15,6 @@
 import argparse
 import re
 import logging
-from readpaf import parse_paf
 from sys import stdin
 from pyfaidx import Fasta
 from pyfastx import Fastq

src/realign.py

@@ -1,6 +1,6 @@
 import os
 import argparse
-from readpaf import parse_paf
+from src import readpaf_local
 import pysam
 
 BAM_CMATCH, BAM_CINS, BAM_CDEL, BAM_CREF_SKIP, BAM_CSOFT_CLIP, BAM_CHARD_CLIP, BAM_CPAD, BAM_CEQUAL, BAM_CDIFF, BAM_CBACK = range(10)
@@ -30,7 +30,7 @@ def run(args):
     # inefficient
     paf_file = open(args.paf, "r")
     paf_dic = {}
-    for record in parse_paf(paf_file):
+    for record in readpaf_local.parse_paf(paf_file):
         paf_dic[record.query_name] = record
 
     processed_sam_record_count = 0

src/reform.py

@@ -1,6 +1,5 @@
 import os
 import argparse
-from readpaf import parse_paf
 import pysam
 
 DEFAULT_KMER_SIZE = 9