Merge pull request #5 from gallardoalba/add_hifiasm_bloom_filter_para…

…meter

Include HiFiasm bloom filter option

README.md

@@ -107,7 +107,7 @@ This command will output you NC_016067.1.fasta and NC_016067.1.gb that you will
 ```
 Usage: 'python mitohifi_v2.py -r "f1.fasta f2.fasta f3.fasta" -f reference.fasta -g reference.gb  -t <int> -o <int> '
 
-usage: mitohifi_v2.py (-r R | -c C) [-h] -f F -g G -t T [-p P]
+usage: mitohifi_v2.py (-r R | -c C) [-h] -f F -g G -t T [-p P] [-m M]
                       [--circular-size CIRCULAR_SIZE]
                       [--circular-offset CIRCULAR_OFFSET] [-o O]
 
@@ -123,6 +123,7 @@ Arguments:
                         hifiams, minimap2, samtools and blast
   -p P                  -p: Percentage of query in the blast match with close-
                         related mito
+  -m M                  -m: Number of bits for HiFiasm bloom filter [it maps to -f in HiFiasm] (default = 0)
   --circular-size CIRCULAR_SIZE
                         Size to consider when checking for circularization
   --circular-offset CIRCULAR_OFFSET

getReprContig.py

@@ -42,12 +42,11 @@ def get_largest_cluster(cdhit_clstr_file):
                 #clusters[cluster_id] = 0
             else:
                 curr_sequences.append(line.strip())
-        # catch the last cluster        
+        # catch the last cluster
         if len(curr_sequences) > largest_cluster_len:
             largest_cluster = curr_cluster
             largest_cluster_len = len(curr_sequences)
             largest_cluster_seqs = curr_sequences
-
     for sequence in largest_cluster_seqs:
         if sequence[-1] == "*":
             representative_seq = sequence

mitohifi_v2.py

@@ -106,6 +106,7 @@ def main():
     parser.add_argument("-g", help= "-k: Close-related species Mitogenome in genebank format", required = "True")
     parser.add_argument("-t", help= "-t: Number of threads for (i) hifiasm and (ii) the blast search", required = "True", type=int)
     parser.add_argument("-p", help="-p: Percentage of query in the blast match with close-related mito", type=int, default=50)
+    parser.add_argument("-m", help="-m: Number of bits for HiFiasm bloom filter [it maps to -f in HiFiasm] (default = 0)", type=int, default=0)
     parser.add_argument('--circular-size', help='Size to consider when checking for circularization', type=int, default=220)
     parser.add_argument('--circular-offset', help='Offset from start and finish to consider when looking for circularization', type=int, default=40)
     parser.add_argument("-o", help="""-o: Organism genetic code following NCBI table (for mitogenome annotation):
@@ -151,7 +152,7 @@ def main():
         print("\nNow let's run hifiasm to assemble the mapped and filtered reads!\n")
 
         with open("hifiasm.log", "w") as hifiasm_log_f:
-            subprocess.run(["hifiasm", "-t", str(args.t), "-o", "gbk.HiFiMapped.bam.filtered.assembled", "gbk.HiFiMapped.bam.filtered.fasta", ], stderr=subprocess.STDOUT, stdout=hifiasm_log_f)
+            subprocess.run(["hifiasm", "-t", str(args.t), "-f", str(args.m), "-o", "gbk.HiFiMapped.bam.filtered.assembled", "gbk.HiFiMapped.bam.filtered.fasta", ], stderr=subprocess.STDOUT, stdout=hifiasm_log_f)
 
         gfa2fa_script = os.path.join(os.path.dirname(__file__),"gfa2fa") # gets path to gfa2fa script
 
@@ -173,7 +174,7 @@ def main():
         fixContigHeaders.fix_headers(original_contigs, "fixed_header_contigs.fasta")
 
         os.remove(original_contigs) # remove original contig file  
-        os.rename("fixed_header_contigs.fasta", original_contigs) # replace original contigs file by the version that has the headers fixed
+        shutil.move("fixed_header_contigs.fasta", original_contigs) # replace original contigs file by the version that has the headers fixed
 
         contigs = original_contigs