software to identify primers that can be used to distinguish genomes
If you use this software, please cite our article in the Journal of Open Source Software.
Note
primerForge is incompatible with versions python3.8 and below and python3.12 and above.
primerForge can be installed with the following commands:
pip install primerforge
conda install ispcrNote
Pixi is experimental
To install with Pixi, download pixi.toml and then run the following commands:
pixi shellA conda installation is currently unavailable due to khmer being unsupported. We are actively working to resolve this.
Note
This might take up to ten minutes.
git clone https://github.com/dr-joe-wirth/primerForge.git
conda env create -f primerForge/environment.yml
conda activate primerforgeA Docker image for the latest release is available at DockerHub
If primerForge is installed correctly, then the following command should execute without errors:
primerForge --check_installIf you installed manually, you may need to use the following command instead:
python primerForge.py --check_installIn order to run unit tests, install primerForge using the instructions above. You will also need to clone the repository if you haven't already:
git clone https://github.com/dr-joe-wirth/primerForge.gitOnce installed and cloned, run the following commands to run the unit tests:
python3 -m unittest discover -s ./primerForge/bin/unit_tests/ -p "*_test.py"usage:
primerForge [-ioBubfpgtrdnkvh]
required arguments:
-i, --ingroup [file] ingroup filename or a file pattern inside double-quotes (eg."*.gbff")
optional arguments:
-o, --out [file] output filename for primer pair data (default: results.tsv)
-B, --bed_file [file] output filename for primer data in BED file format (default: primers.bed)
-u, --outgroup [file] outgroup filename or a file pattern inside double-quotes (eg."*.gbff")
-b, --bad_sizes [int,int] a range of PCR product lengths that the outgroup cannot produce (default: same as '--pcr_prod')
-f, --format [str] file format of the ingroup and outgroup [genbank|fasta] (default: genbank)
-p, --primer_len [int(s)] a single primer length or a range specified as 'min,max'; (minimum 10) (default: 16,20)
-g, --gc_range [float,float] a min and max percent GC specified as a comma separated list (default: 40.0,60.0)
-t, --tm_range [float,float] a min and max melting temp (°C) specified as a comma separated list (default: 55.0,68.0)
-r, --pcr_prod [int(s)] a single PCR product length or a range specified as 'min,max' (default: 120,2400)
-d, --tm_diff [float] the maximum allowable Tm difference °C between a pair of primers (default: 5.0)
-n, --num_threads [int] the number of threads for parallel processing (default: 1)
-k, --keep keep intermediate files (default: False)
-v, --version print the version
-h, --help print this message
--check_install check installation
--debug run in debug mode
--advanced print advanced options
--primer3_mv_conc [float] monovalent cation concentration (mM) (default: 50.0)
--primer3_dv_conc [float] divalent cation concentration (mM) (default: 1.5)
--primer3_dntp_conc [float] dNTP concentration (mM) (default: 0.6)
--primer3_dna_conc [float] template DNA concentration (nM) (default: 50.0)
--primer3_temp_c [float] simulation temp (°C) for ΔG (default: 37.0)
--primer3_max_loop [int] maximum size (bp) of loops in primer secondary structures (default: 30)
--isPcr_minGood [int] minimum size (bp) where there must be 2 matches for each mismatch (default: 6)
--isPcr_minPerfect [int] minimum size (bp) of perfect match at 3' end of primer (default: 8)
--isPcr_tileSize [int] the size of match that triggers an alignment (default: 10)
--temp_tolerance [float] minimum number of degrees (°C) below primer Tm allowed for secondary structure Tm (default: 5.0)
--max_repeats [int] maximum allowed length (bp) of homopolymers (repeats) in primer sequences (default: 4)
--bin_size [int] maximum allowed length (bp) of contiguous regions of overlapping primers (bins) (default: 64)
The results produced by primerForge list one primer pair per line, and these pairs are sorted based on the following criteria (in order of importance):
- largest difference from the ingroup PCR product range for outgroup PCR products
- lowest number of outgroup PCR products
- lowest G+C difference between the pair
- lowest Tm difference between the pair
- lowest Tm for heterodimers
- lowest Tm (sum) for homodimers
- lowest Tm (sum) for hairpins
- lowest variance in ingroup PCR product sizes
- largest median ingroup PCR product size
See the example for details about the file format.
The BED file produced by primerForge has the following format:
| column number | title | definition |
|---|---|---|
| 1 | contig | the name of the contig |
| 2 | start | the start coordinate of the primer |
| 3 | end | the end coordinate of the primer (exclusive) |
| 4 | sequence | the sequence of the primer |
| 5 | pair number | an integer that links forward and reverse primers; overloading the "score" field traditionally found in BED file format |
| 6 | strand | the DNA strand ('+' or '-') |
flowchart TB
ingroup[/"ingroup genomes"/]
ingroup --> A
%% get unique kmers
subgraph A["for each genome"]
uniqKmer["get unique kmers"]
end
%% get shared kmers
sharedKmers(["shared kmers"])
uniqKmer -- intersection --> sharedKmers
%% get candidate kmers
subgraph B["for each genome"]
subgraph B0["for each kmer start position"]
subgraph B1["pick one kmer"]
GC{"GC in
range?"}
Tm{"Tm in
range?"}
homo{"repeats
≤ 3bp?"}
hair{"no hairpins?"}
dime{"no homo-
dimers?"}
GC-->Tm-->homo-->hair-->dime
end
end
end
%% connections up to candidate kmers
sharedKmers --> B
dump1[/"checkpoint"/]
sharedKmers --> dump1
candidates(["unique, shared kmers; one per start position"])
dime --> candidates
%% get primer pairs
subgraph C["for one genome"]
bin1["bin overlapping kmers (64bp max)"]
bin2["remove kmers that are
substrings of other kmers"]
bin3["get bin pairs"]
%% evaluate one kmer pair
subgraph C0["for each bin pair"]
size{"is PCR
size ok?"}
subgraph C1["for each primer pair"]
prime{"is 3' end
G or C?"}
temp{"is Tm
difference ok?"}
gc2{"is GC
difference ok?"}
hetero{"no hetero-
dimers?"}
end
size --> C1
end
end
candPair(["candidate primer pairs"])
allSharePair(["all shared primer pairs"])
%% get shared primer pairs
subgraph D["for each candidate primer pair"]
subgraph D0["for each other genome"]
pcr{"is PCR
size ok?"}
end
end
bin1 --> bin2
bin2 --> bin3
bin3 --> C0
prime --> temp --> gc2 --> hetero --> candPair
candPair --> D
pcr --> allSharePair
%% one pair per bin pair
subgraph E["for each bin pair"]
keep["keep only one primer pair"]
end
selectedSharePair(["selected shared
primer pairs"])
dump2[/"checkpoint"/]
dump3[/"checkpoint"/]
candidates --> dump2
candidates --> C
allSharePair --> E
keep --> selectedSharePair
selectedSharePair --> dump3
%% outgroup removal
outgroup[/"outgroup genomes"/]
subgraph F["for each outgroup genome"]
subgraph F0["for each primer pair"]
ogsize{"PCR size outside
disallowed range?"}
end
end
selectedSharePair --> F0
outgroup --> F
noout(["primer pairs absent from outgroup"])
ogsize --> noout
noout --> dump4
dump4[/"checkpoint"/]
ispcr["filter pairs using isPcr"]
final(["final primer pairs"])
dump5[/"checkpoint"/]
write[/"sort pairs and
write to files"/]
noout --> ispcr
ispcr --> final
final --> dump5
final --> write
This example assumes you have already installed primerForge as described above.
In this example, we will use primerForge to find pairs of primers between 18bp and 24bp that can be used to differentiate three strains of Mycoplasma mycoides subspecies mycoides (the "ingroup") from two strains of Mycoplasma mycoides subspecies capri (the "outgroup"). The primer pairs identified by primerForge are predicted to produce a single PCR product between 500bp and 2000bp in the ingroup. These same primer pairs are predicted to produce PCR products <300bp, >3000bp, or no PCR products at all in the outgroup.
In order to get started, create a directory called mycoplasma_test and move into it:
mkdir ./mycoplasma_test
cd ./mycoplasma_testNext, download the following Mycoplasma mycoides genomes using the following commands:
curl ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/034/305/GCF_003034305.1_ASM303430v1/GCF_003034305.1_ASM303430v1_genomic.gbff.gz | gzip -d > ./i1.gbff
curl ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/034/275/GCF_003034275.1_ASM303427v1/GCF_003034275.1_ASM303427v1_genomic.gbff.gz | gzip -d > ./i2.gbff
curl ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/034/345/GCF_003034345.1_ASM303434v1/GCF_003034345.1_ASM303434v1_genomic.gbff.gz | gzip -d > ./i3.gbff
curl ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/900/489/555/GCF_900489555.1_MMC68/GCF_900489555.1_MMC68_genomic.gbff.gz | gzip -d > ./o1.gbff
curl ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/018/389/745/GCF_018389745.1_ASM1838974v1/GCF_018389745.1_ASM1838974v1_genomic.gbff.gz | gzip -d > ./o2.gbffIf you cannot download the genbank files using curl, you can download them manually from NCBI by replacing ftp:// with http:// and copying and pasting each address into your web browser (eg. http://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/034/305/GCF_003034305.1_ASM303430v1/GCF_003034305.1_ASM303430v1_genomic.gbff.gz) and then using gzip -d on the downloaded file to uncompress it. Finally, be sure to rename each file as shown above (eg. mv GCF_003034305.1_ASM303430v1_genomic.gbff i1.gbff).
We will use the following flags to specify specific parameters for this example:
- The
--ingroupand--outgroupflags are both file patterns for the ingroup and outgroup genomes, respectively. It is important that this pattern is enclosed in double-quotes as shown below. - The
--pcr_prodflag indicates what sizes we want for ingroup products (500bp to 2,000bp) - The
--bad_sizesflag indicates what sizes we do not want for outgroup products (300bp to 3,000bp). - The
--primer_lenflag indicates what length our primers can be (18bp to 24bp)
You can get a list of all available flags using the command primerForge --help.
Run primerForge using the following command (requires at least 2 Gb of RAM):
primerForge --ingroup "./i*gbff" --outgroup "./o*gbff" --pcr_prod 500,2000 --bad_sizes 300,3000 --primer_len 18,24After running the command, you should see something like this printed to the screen:
dumping Parameters to '_pickles/parameters.p' ... done 00:00:00.49
identifying kmers suitable for use as primers in all 3 ingroup genome sequences
getting shared ingroup kmers that appear once in each genome ... done 00:01:36.15
dumping shared kmers to '_pickles/sharedKmers.p' ... done 00:00:11.22
evaluating 2430140 kmers ... done 00:01:51.73
identified 30413 candidate kmers
dumping candidate kmers to '_pickles/candidates.p' ... done 00:00:01.55
done 00:03:42.94
identifying pairs of primers found in all ingroup sequences ... done 00:00:11.62
identified 16050 primer pairs shared in all ingroup sequences
dumping unfiltered pairs to '_pickles/pairs.p' ... done 00:00:00.56
removing primer pairs present in the outgroup sequences
getting outgroup PCR products ... done 00:00:01.03
filtering primer pairs ... done 00:00:00.54
processing outgroup results ... done 00:00:00.54
removed 5905 pairs present in the outgroup (10145 pairs remaining)
dumping filtered pairs to '_pickles/pairs_noOutgroup.p' ... done 00:00:00.54
validating primer pairs with isPcr ... done 00:00:03.44
removed 3659 pairs not validated by isPcr (6486 pairs remaining)
dumping validated pairs to '_pickles/pairs_noOutgroup_validated.p' ... done 00:00:00.54
writing 6486 primer pairs to 'results.tsv' ... done 00:00:11.62
total runtime: 00:04:14.50
As we can see, primerForge found 30,413 kmers that were suitable for use as a primer in all three ingroup genomes. It then went on to identify 16,050 primer pairs that would produce PCR products between 500bp and 2000bp in the ingroup genomes. Next, it found that of those 16,050 pairs, 5,905 of them formed PCR products between 300bp and 3000bp in one or more of the outgroup genomes. Finally, it used isPcr to validate the remaining 10,145 primer pairs resulting in 6,486 primer pairs being written to file.
primerForge generated results.tsv, the file that contains the sequences and details for each primer pair, and primerForge.log. Here are a few lines from results.tsv:
| fwd_seq | fwd_Tm | fwd_GC | rev_seq | rev_Tm | rev_GC | i1.gbff_contig | i1.gbff_length | i2.gbff_contig | i2.gbff_length | i3.gbff_contig | i3.gbff_length | o1.gbff_contig | o1.gbff_length | o2.gbff_contig | o2.gbff_length |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| AGAAGCAAAGGATATGGGAACAAC | 57.1 | 41.7 | AAATCAACACCCTCAATAAGCTCC | 57.1 | 41.7 | NZ_LAUX01000078.1 | 804 | NZ_LAUV01000035.1 | 804 | NZ_LAUY01000092.1 | 804 | NA | 0 | NA | 0 |
| TCCATCTAATGAAGATCAACCAGG | 55.9 | 41.7 | CCCTAATTGTGATGAGTTGACAAC | 55.9 | 41.7 | NZ_LAUX01000010.1 | 710 | NZ_LAUV01000042.1 | 713 | NZ_LAUY01000078.1 | 710 | NA | 0 | NA | 0 |
| CATCAGCTGTTGTAAATAACCCAC | 56.2 | 41.7 | GTGGAGCTATGAAACCAATATCAG | 55.3 | 41.7 | NZ_LAUX01000117.1 | 1694 | NZ_LAUV01000064.1 | 1694 | NZ_LAUY01000018.1 | 1694 | NZ_LS483503.1 | 3212 | NA | 0 |
The first six columns show the forward and reverse sequences (5' --> 3') as well as their melting temperatures and their G+C content (mol %). Next, for each genome it lists the contig and the PCR product size that is predicted be produced by this pair. For example, the first pair of primers are predicted to produce PCR products of 804bp the ingroup genomes, and the binding sites for this primer pair in the files i1.gbff, i2.gbff, and i3.gbff can be found in contigs NZ_LAUX01000078.1, NZ_LAUV01000035.1, and NZ_LAUY01000092.1, respectively. This same pair is not predicted to produce any PCR products in either outgroup genome. Similarly, the third pair is predicted to produce a PCR product size of 1,694bp in all three ingroup genomes, no products o2.gbff, and a PCR product >3,000bp in o1.gbff.
If a primer pair is predicted to produce multiple products in an outgroup genome, then the contig column and the size column will list contigs and sizes in a comma-separated list linked by position. For example, if a primer pair was expected to produce a product of 1,990bp in contig_1 and 2,024bp in contig_2 in the genome file example.gbff, then the columns for this genome would look like this:
| example.gbff_contig | example.gbff_length |
|---|---|
| 1990,2024 | contig1,contig2 |
Note
Multiple PCR products will only ever be predicted for outgroup genomes as primerForge does not allow such primer pairs in the ingroup genome.
This error occurs if you specify a wildcard representing input files without enclosing them in quotes. For example, this will cause the error:
primerForge --ingroup ./i*gbffand this will fix it:
primerForge --ingroup "./i*gbff"The same holds true for the --outgroup flag.
This error can also occur if you have inadvertently included a space in any of the arguments passed to other flags.
This error occurs if the specified file(s) cannot be found or if the file format does not match the --format flag (default = genbank). Check that the file path is correct and the files can be read. If they are correct, then double check that you have specified --format fasta if you are working with fasta files.
This error occurs if primerForge cannot find kmers that are shared in all the ingroup genomes. This can occur if the input genomes are too distantly related, or if one or more of the genomes is of very poor quality. To diagnose this, try repeating the command but include the --debug flag. This will report which ingroup genome is causing the number of shared kmers to drop to zero in the primerForge.log file.
This error occurs if all the primer pairs primerForge identified cannot be used to distinguish the ingroup from the outgroup. This most often occurs because one or more members of the outgroup is too closely-related to the ingroup. To diagnose this, try repreating the command but include the --debug flag. This will report which outgroup genome is causing the number of shared kmers to drop to zero in the primerForge.log file. Alternatively, you can expand your search by widening the ranges passed to the flags --pcr_prod and/or --bad_sizes.
Thank you for your interest in contributing! Please see the contributing guidelines for more information.