add getpos_primers.py to add positions for primerdesign.py output

Author: Runsheng

bin/getpos_primers.py

@@ -0,0 +1,83 @@
+#!/usr/bin/env python
+#-*- coding: utf-8 -*-
+#@Time    : 17/7/2023 10:18 pm
+#@Author  : Runsheng
+#@File    : getpos_primers.py
+
+"""
+get the position using ispcr function for all output from primerdiffer.py
+run example:
+cd /t1/ref_BN
+getpos_primers.py -f primer10.txt -g cb5.fa -o primer10_pos.txt
+"""
+
+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", "--file",
+                    help="the 4 col table generated by primerdesign.py")
+
+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="primerpos.txt",
+                    help="output file, contains 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
+seqdic=dic2dic(fasta2dic(args.genome))
+
+def read_primerdiffer_out(filename):
+    lines=[]
+    with open(filename, "r") as f:
+        for line in f.readlines():
+            line_l=line.strip().split("\t")
+            lines.append(line_l)
+    return lines
+
+lines=read_primerdiffer_out(args.file)
+
+lines_new=[]
+for line_l in lines:
+    name_raw, forward, reverse, length_str=line_l
+
+    product_l=insilicon_pcr(primer_left=forward,
+              primer_right=reverse,
+              db=args.genome,
+              cutoff_alignlength=args.alignlen,
+              cutoff_free3=args.free3len,
+              product_cutoff=args.productlen
+              )
+
+    outdic=product2seqdic(product_l,seqdic)
+    name_pos=list(outdic.keys())[0]
+    seq=outdic[name_pos]
+    line_one=[name_raw, forward, reverse, length_str, name_pos, seq]
+    lines_new.append(line_one)
+
+with open(args.out, "w") as fw:
+    for line_one in lines_new:
+        fw.write("\t".join(line_one))
+        fw.write("\n")

primerdiffer/__init__.py

@@ -1 +1 @@
-__version__= "0.1.5"
+__version__= "0.1.6"

readme.md

@@ -70,16 +70,47 @@ primerdesign.py -g1 cb5.fa -g2 cn3_new.fa -pos "ChrX:12881200-15106660" --interv
 ```
 
 ```
-# check the result in file "primers_ChrX:12881200-15106660.txt"
-head primers_ChrX\:12881200-15106660.txt
+# check the result in file "primers_ChrX_12881200-15106660.txt"
+head primers_ChrX_12881200-15106660.txt
 #ChrX:12881200-12881700	GATCCAAAACATGAGTGGCC	CGAGATCATTGGCTCAAAGT	287
 #ChrX:12886200-12886700	GTTTTCTCTTCAAGTGCCCG	CTCCCACATCTTGTAGGTCC	416
 #ChrX:12891200-12891700	GTAGATGCTGTTGAGGCTCT	CGAGTGGGACATTGTCAGTA	299
 #ChrX:12896200-12896700	GGCGCATTATACGAAGCTTT	TTCCTGCTGCCAGATAGAAG	362
 ```
 
+## Case example: use getpos_primers.py to annotate the primerdiffer.py output with position and in silicon PCR result
 
-Use in silico PCR to get the position and the product of the primer
+```bash
+usage: getpos_primers.py [-h] [-d WKDIR] [-f FILE] [-g GENOME] [--alignlen ALIGNLEN] [--free3len FREE3LEN] [--productlen PRODUCTLEN] [-o OUT]
+
+optional arguments:
+  -h, --help            show this help message and exit
+  -d WKDIR, --wkdir WKDIR
+                        The dir path contain the file, if not given, use the current dir
+  -f FILE, --file FILE  the 4 col table generated by primerdesign.py
+  -g GENOME, --genome GENOME
+                        the fasta file used to design primer
+  --alignlen ALIGNLEN   the cutoff of primer min align length as a right hit, default is 16
+  --free3len FREE3LEN   the cutoff of primer 3' align length as a right hit, default is 2
+  --productlen PRODUCTLEN  the cutoff of max product which will be treated as a false priming, default is 2000.
+  -o OUT, --out OUT     output file, contains possible amplification regions of this primer pair
+```
+
+```bash
+# example:  first check the input file format, the col4 is from the primerdesign.py output
+head primers_ChrX_12881200-15106660.txt 
+#ChrX:12881200-12881700	GATCCAAAACATGAGTGGCC	CGAGATCATTGGCTCAAAGT	287
+
+# run get_posprimers.py
+getpos_primers.py -f primers_ChrX_12881200-15106660.txt -g cb5.fa -o primer_pos.txt
+
+# the output is 6 col table
+#ChrX:12881200-12881700	GATCCAAAACATGAGTGGCC	CGAGATCATTGGCTCAAAGT	287 ChrX:12881301-12881587  ATCCAAAACATGAG....
+```
+
+
+## Case example: in silico PCR for one primer
+Use in silico PCR to get the position and the product of one primer
 ```
 usage: ispcr.py [-h] [-d WKDIR] [-f FORWARD] [-r REVERSE] [-g GENOME] [--alignlen ALIGNLEN]
                 [--free3len FREE3LEN] [--productlen PRODUCTLEN] [-o OUT]
@@ -102,7 +133,6 @@ optional arguments:
   -o OUT, --out OUT     output file, contains all possible amplification regions of this primer pair
 ```
 
-## Case example: in silico PCR product
 
 ```bash
 ispcr.py -f gcactttcatgtccctcaac -r cactctattctcaccccacc -g cb5.fa -o ispcr.fa
@@ -119,6 +149,7 @@ head ispcr.fa
 #ATGATGATGATGATGGTGGGGTGAGAATAGAGT
 ```
 
+
 ## Case example: Design general purpose primers with given parameters
 The default primer design parameter is described in [general_setting.py](https://github.com/Runsheng/primerdiffer/blob/master/primerdiffer/general_settings.py):
 ```python
@@ -163,5 +194,4 @@ Please kindly cite our paper in bioinformatics if you find primerdiffer or prime
 
 ## Other functions:
 1. To design the primer using VCF file for closely related haplotypes (i.e., strain/individual level differences). Please check primervcf package.
-2. Todo: To update the RFLP method for primer design to differ sequences with almost identical sequence.
-
+2. add the primerdffer output table description and ispcr function for the output.

setup.py

@@ -18,7 +18,8 @@
     install_requires = ["primer3-py>=0.6.1",
                         "biopython>=1.78"],
     scripts = ['bin/primerdesign.py',
-               'bin/ispcr.py'
+               'bin/ispcr.py',
+               'bin/getpos_primers.py'
                 ]
 )
 

test/test_walk_chr.py

@@ -30,8 +30,8 @@ def setUp(self):
     def test_config_dump_load(self):
         import json
         from primerdiffer.general_settings import primer3_general_settings
-        with open("./test/primer3config", "w") as fw:
-            json.dump(primer3_general_settings, )
+        with open("./primer3config", "w") as fw:
+            json.dump(primer3_general_settings, fw)
 
 
     def test_walk_chr_dense(self):