add scripts and finalize v0.1.1

Author: Runsheng

bin/__init__.py

@@ -0,0 +1,64 @@
+#!/usr/bin/env python
+#-*- coding: utf-8 -*-
+#@Time    : 9/10/2022 5:25 PM
+#@Author  : Runsheng
+#@File    : ispcr.py
+
+"""
+in silico PCR script
+"""
+
+import argparse
+import os,sys
+
+# self import
+from primerdiffer.primer_check import insilicon_pcr, product2seqdic
+from primerdiffer.utils import checkblastdb, dic2fasta, dic2dic,fasta2dic
+
+
+parser=argparse.ArgumentParser()
+parser.add_argument("-d", "--wkdir", default=None,
+                    help="The dir path contain the file, if not given, use the current dir")
+# input
+parser.add_argument("-f", "--forward",
+                    help="the forward primer sequence")
+parser.add_argument("-r", "--reverse",
+                    help="the reverse primer sequence")
+
+parser.add_argument("-g", "--genome",
+                    help="the fasta file used to design primer")
+
+# primer cutoff
+parser.add_argument("--alignlen", type=int, default=16,
+                    help="the cutoff of primer min align length as a right hit, default is 16")
+parser.add_argument("--free3len",  type=int, default=2,
+                    help="the cutoff of primer 3' align length as a right hit, default is 2")
+parser.add_argument("--productlen",  type=int, default=2000,
+                    help="the cutoff of max product which will be treated as a false priming, default is 2000.")
+
+# output paramenter
+parser.add_argument("-o", "--out", default="product.fasta",
+                    help="output file, contains all possible amplification regions of this primer pair")
+
+# default handler
+args = parser.parse_args(args=None if sys.argv[1:] else ['--help'])
+wkdir=os.getcwd() if args.wkdir is None else args.wkdir
+
+# main
+checkblastdb(args.genome) # check blastdb, if no, run makeblastdb
+
+product_l=insilicon_pcr(primer_left=args.forward,
+              primer_right=args.reverse,
+              db=args.genome,
+              cutoff_alignlength=args.alignlen,
+              cutoff_free3=args.free3len,
+              product_cutoff=args.productlen
+              )
+seqdic=dic2dic(fasta2dic(args.genome))
+outdic=product2seqdic(product_l,seqdic)
+dic2fasta(outdic, args.out)
+
+
+
+
+

bin/ispcr.py

@@ -0,0 +1,72 @@
+#!/usr/bin/env python
+#-*- coding: utf-8 -*-
+#@Time    : 9/10/2022 4:30 PM
+#@Author  : Runsheng
+#@File    : primerdesign.py
+
+"""
+The main script used to run cmd orders
+"""
+
+import argparse
+import os,sys
+
+#currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
+#parentdir = os.path.dirname(os.path.dirname(currentdir))
+#sys.path.insert(0,parentdir)
+
+# self import
+from primerdiffer.walk_chr import flow_walk_chr
+
+parser=argparse.ArgumentParser()
+parser.add_argument("-d", "--wkdir", default=None,
+                    help="The dir path contain the file, if not given, use the current dir")
+# input
+parser.add_argument("-g1", "--genome1",
+                    help="the fasta file used to design primer")
+parser.add_argument("-g2", "--genome2",
+                    help="the fasta file used to check false priming")
+parser.add_argument("-pos", "--position",default=None,
+                    help="position on genome1 to design primer, use string with IGV format, like ChrX:1956230-1976220")
+
+# primer cutoff
+parser.add_argument("--alignlen", type=int, default=16,
+                    help="the cutoff of primer min align length as a right hit, default is 16")
+parser.add_argument("--free3len",  type=int, default=2,
+                    help="the cutoff of primer 3' align length as a right hit, default is 2")
+parser.add_argument("--productlen",  type=int, default=2000,
+                    help="the cutoff of max product which will be treated as a false priming, default is 2000.")
+
+# hit cutoff
+parser.add_argument("-h1","--hit1", type=int, default=1,
+                    help="the cutoff of max number of in-silicon PCR product which can be found in genome1. default is 1")
+parser.add_argument("-h2","--hit2",  type=int, default=0,
+                    help="the cutoff of max number of in-silicon PCR product which can be found in genome2, default is 0")
+
+# run parameter, how dense the primer should be
+parser.add_argument("-i","--interval", type=int, default=5000,
+                    help="interval is the region begins to pick primers, default is 5000. If 5k is the unit, will pick one primer each 5k")
+parser.add_argument("-j","--jump",  type=int, default=500,
+                    help="jump is the region to jump inside intervals if the prvious interval can not generate a valid primer"
+                         ", the smaller, more sites to check. Default is 500.")
+# output paramenter
+parser.add_argument("--prefix", default="primers",
+                    help="prefix of output file, default is primers")
+
+
+# default handler
+args = parser.parse_args(args=None if sys.argv[1:] else ['--help'])
+wkdir=os.getcwd() if args.wkdir is None else args.wkdir
+
+flow_walk_chr(wkdir=wkdir,
+              genome1=args.genome1,
+              genome2=args.genome2,
+              pos_str=args.position,
+              cutoff_alignlength=args.alignlen,
+              cutoff_free3=args.free3len,
+              product_cutoff=args.productlen,
+              db1_maxhit=args.hit1,
+              db2_maxhit=args.hit2,
+              interval=args.interval,
+              jump=args.jump,
+              out_prefix=args.prefix)

bin/primerdesign.py

@@ -1 +1 @@
-__version__= "0.1.0"
+__version__= "0.1.1"

primerdiffer/__init__.py

@@ -4,18 +4,20 @@
 # @Author  : Runsheng
 # @Email   : Runsheng.lee@gmail.com
 # @File    : primer_check.py
-import primer3
 
-from primerdiffer.general_settings import primer3_general_settings
 
 try:
     from StringIO import StringIO ## for Python 2
 except ImportError:
     from io import StringIO ## for Python 3
 
+import primer3
 from Bio.Blast.Applications import NcbiblastnCommandline
 from Bio.Blast import NCBIXML
 
+from primerdiffer.general_settings import primer3_general_settings
+from primerdiffer.utils import fasta2dic,dic2dic,chr_select, reverse_complement
+
 
 def my_design_primer(name,seq,primer3_settings=primer3_general_settings):
     """
@@ -66,7 +68,7 @@ def filter_hsp(blast_records,query,cutoff_alignlength=16,cutoff_free3=2, debugmo
     return keep
 
 
-def insilicon_pcr(primer_left, primer_right, db, cutoff_alignlength=16, cutoff_free3=2, profuct_cutoff=2000,
+def insilicon_pcr(primer_left, primer_right, db, cutoff_alignlength=16, cutoff_free3=2, product_cutoff=2000,
                   debugmod=False):
     """
     para: the left and right primers
@@ -85,12 +87,29 @@ def insilicon_pcr(primer_left, primer_right, db, cutoff_alignlength=16, cutoff_f
         for pr in p_right:
             # print pl, pr
             # use only the start to get a approx length, also, the direction for left and right primer should be different
-            if pl[0] == pr[0] and abs(pl[1] - pr[1]) <= profuct_cutoff and pl[-1] * pr[-1] == -1:
+            if pl[0] == pr[0] and abs(pl[1] - pr[1]) <= product_cutoff and pl[-1] * pr[-1] == -1:
                 possible_product.append((pl[0], pl[1], pr[1]))
 
     return possible_product
 
 
+def product2seqdic(product_l, seqdic):
+    """
+    use one line from product_l to write the sequence
+    """
+    outdic={}
+    for product_1 in product_l:
+        chro, start, end=product_1
+        # the product could be forward (as usual) or reverse (F as R and R as F)
+        if end>=start:
+            name, seq=chr_select(seqdic, chro, start, end)
+            outdic[name]=seq
+        else:
+            name, seq=chr_select(seqdic, chro, end, start)
+            outdic[name+"_RC"]=reverse_complement(seq)
+    return outdic
+
+
 def _is_nofalse_primer(blast_records,query,debugmod=False):
     """
     :param blast_records: input a blast record in XML format

primerdiffer/primer_check.py

@@ -10,8 +10,21 @@
 """
 
 # Third oart import
-from Bio import SeqIO
+from os.path import exists
 
+from Bio import SeqIO, Seq
+from Bio.Blast.Applications import NcbimakeblastdbCommandline
+
+
+def reverse_complement(seq):
+    """
+        Given: A DNA string s of length at most 1000 bp.
+        Return: The reverse complement sc of s.
+        due to the complement_map,
+        the symbol such as \n and something else is illegal
+        the input need to be pure sequence
+    """
+    return str(Seq.reverse_complement(seq))
 
 def fasta2dic(fastafile):
     """
@@ -35,6 +48,13 @@ def dic2dic(record_dict):
         seq_dict[k]=seq.upper()
     return seq_dict
 
+def dic2fasta(record_dict, out):
+    with open(out, "w") as fw:
+        for k, v in record_dict.items():
+            fw.write(">"+k+"\n")
+            fw.write(v)
+            fw.write("\n")
+
 
 def tuple_to_pos_str(pos_t):
     """
@@ -68,3 +88,21 @@ def chr_select(seq_dict, chro, start,end):
     return name,seq
 
 
+def makeblastdb(genomefile):
+    cline = NcbimakeblastdbCommandline(dbtype="nucl",
+                                       input_file=genomefile)
+    NcbimakeblastdbCommandline(cmd='makeblastdb', dbtype='prot', input_file='NC_005816.faa')
+    print(cline)
+    cline()
+
+
+def checkblastdb(genomefile):
+    """
+    check if the blastdb exist, if not, create one
+    """
+    dbfile=genomefile+".nsq"
+    if exists(dbfile):
+        return 0
+    else:
+        makeblastdb(genomefile)
+        return 1

primerdiffer/utils.py

@@ -9,16 +9,15 @@
 
 # self import
 import os
-from os.path import exists
-
-from Bio.Blast.Applications import NcbimakeblastdbCommandline
 
 from primerdiffer.primer_check import my_design_primer, primer_check
-from primerdiffer.utils import dic2dic, fasta2dic, chr_select, tuple_to_pos_str, pos_str_to_tuple
+from primerdiffer.utils import dic2dic, fasta2dic, chr_select, tuple_to_pos_str, pos_str_to_tuple, checkblastdb
 
 
 def walk_chr_dense(genome, chro, start, end, db1, db2,
-                   interval=500000, jump=4000, out_prefix="primers"):
+                   cutoff_alignlength=16, cutoff_free3=2, product_cutoff=2000,
+                   db1_maxhit=1, db2_maxhit=0,
+                   interval=500000, jump=4000, out_prefix="primers", debugmod=False):
     """
     :param genome: genome is a dict in name:seq,
     :param chro:
@@ -43,7 +42,10 @@ def walk_chr_dense(genome, chro, start, end, db1, db2,
                 offset += 1
             else:
                 myprimer = my_design_primer(name=name, seq=seq)
-                primer_used = primer_check(myprimer,db1,db2, debugmod=True)
+                primer_used = primer_check(myprimer,db1,db2,
+                                           cutoff_alignlength, cutoff_free3, product_cutoff,
+                                           db1_maxhit, db2_maxhit,
+                                           debugmod=debugmod)
                 if primer_used == 0:
                     offset += 1
                 else:
@@ -62,27 +64,10 @@ def walk_chr_dense(genome, chro, start, end, db1, db2,
     return primer_dict
 
 
-def makeblastdb(genomefile):
-    cline = NcbimakeblastdbCommandline(dbtype="nucl",
-                                       input_file=genomefile)
-    NcbimakeblastdbCommandline(cmd='makeblastdb', dbtype='prot', input_file='NC_005816.faa')
-    print(cline)
-    cline()
-
-
-def checkblastdb(genomefile):
-    """
-    check if the blastdb exist, if not, create one
-    """
-    dbfile=genomefile+".nsq"
-    if exists(dbfile):
-        return 0
-    else:
-        makeblastdb(genomefile)
-        return 1
-
-
-def flow_walk_chr(wkdir, genome1, genome2, pos_str, interval=4000, jump=400, out_prefix="primers"):
+def flow_walk_chr(wkdir, genome1, genome2, pos_str,
+                  cutoff_alignlength=16, cutoff_free3=2, product_cutoff=2000,
+                  db1_maxhit=1, db2_maxhit=0,
+                  interval=4000, jump=400, out_prefix="primers"):
     """
 
     genome1: the genome fasta file used to design primers
@@ -100,10 +85,14 @@ def flow_walk_chr(wkdir, genome1, genome2, pos_str, interval=4000, jump=400, out
 
     # main
     g1=dic2dic(fasta2dic(genome1))
-    #g2=dic2dic(fasta2dic(genome2))
     chro, start, end = pos_str_to_tuple(pos_str)
 
     primer_dict=walk_chr_dense(genome=g1, chro=chro, start=start, end=end,
                                db1=genome1,db2=genome2,
-                   interval=interval, jump=jump, out_prefix=out_prefix)
+                               cutoff_alignlength=cutoff_alignlength,
+                               cutoff_free3=cutoff_free3,
+                               product_cutoff=product_cutoff,
+                               db1_maxhit=db1_maxhit,
+                               db2_maxhit=db2_maxhit,
+                            interval=interval, jump=jump, out_prefix=out_prefix)
     print(primer_dict)