<p>I've found the problem, and it's way less interesting than I initially thought. Essentially, my "eValue" was way too low. Instead of <code>-evalue 100</code>, I'm now using <code>-evalue 5000</code> (Primer-BLAST actually uses 30,000).</p>

<p>For some primers, <code>-evalue 5000</code> gives ~20,000 results. It's beyond the scope of this question, but I'm now filtering on the query/subject Tm for each hit, perfect match at the query 3'-end, and number of query/subject mismatches. My pairwise 
hit comparisons are now detecting plausible off-target amplicons that Primer-BLAST is detecting.</p>

<p>Essentially, I was relying on NCBI's hit ranking system to determine what hits were relevant, but I needed cast a wider net with blastn, and then implement a rational ranking scheme of my own.</p>

CC BY-SA 4.0
2019-03-28T17:54:20.187Z
2019-03-28T22:49:50.650Z
2019-03-28T22:49:50.650Z
mRotten
<p><strong>Summary:</strong></p>

<p>I'm writing an automated tool for designing Sanger sequencing primers (other tools don't do <em>exactly</em> what I want) that will be incorporated in a production diagnostic laboratory for on-demand NGS variant verification. I'm stuck on a problem wherein off-target BLAST hits found by NCBI's Primer-BLAST aren't showing up in my programmatic BLAST query of the same genome assembly. </p>

<p><strong>Background:</strong></p>

<p>After each iteration of design with primer3, I'm using <code>blastn -task blastn</code> to search a local hg38 (patch 12) database for off-target hits. I have logic that does round-robin pairwise comparisons (i.e. compares every hit to every other hit) to determine if each hit pair could result in an amplicon (on same chromosome, opposite-stranded, perfect match at 3'-ends, plus-strand hit is upstream of minus-strand hit, are separated by &lt;= 1000 bp). If my logic finds a viable off-target amplicon, it triggers a total re-design with different design parameters. If no off-targets are found, those primers are finalized, so it's important that I find any and all off-targets.</p>

<p>To make sure that code for viable off-target amplicons is working correctly (that it's triggering re-design or reporting a design failure when appropriate), I'm using NCBI's browser-based Primer-BLAST to verify the absence of viable off-targets, which also uses genome build hg38 patch 12.</p>

<p><strong>The problem:</strong></p>

<p>Primer-BLAST is finding hits that I can't find with my code. Even when I look at every hit returned from the <code>blastn</code> query, I can't find hits that Primer-BLAST is finding. I would like to know if there is some parameter that I am missing, or parameters I could use in addition to those I'm already using, to produce the hits I am getting from Primer-BLAST.</p>

<p>I believe my blast database an my query fasta file are built correctly, because I am getting logical output from the query below.</p>

<p><strong>Programmatic BLAST Query Details:</strong></p>

<p>Query file (query.fasta):</p>

<pre><code>&gt;for_primer
GGATCTGCCTCCTTGCC
&gt;rev_primer
CCCAACAGGGCCCTC
</code></pre>

<p>BLASTn Query Statement:</p>

<pre><code>blastn -db blastn -db /path/to/hg38.fasta -evalue 100 -num_threads 8 -outfmt 6 std gaps nident -dust no -word_size 7 -query /path/to/query.fasta
</code></pre>

<p>The BLAST database used hg38.fa.gz (which I believe to be patch 12) from <a href="http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/" rel="nofollow noreferrer">here</a>. That was unzipped and processed with <code>makeblastdb</code>.</p>

<p><strong>Primer-BLAST Results:</strong></p>

<p>The Primer-BLAST results are <a href="https://www.ncbi.nlm.nih.gov/tools/primer-blast/primertool.cgi?ctg_time=1552618788&amp;job_key=QEqf-GL6b1JIbHVpeAlRWwISQGkvAVt0Lg" rel="nofollow noreferrer">here</a>. To re-create, use the template DNA below, search database = "Genomes for selected organisms...", organism = <em>Homo sapiens</em>, max target size = 1000. The template DNA is:</p>

<pre><code> GAGGAGTCAGGCCTGATCTGAAGGCAGTGGGGTACCCTCTAGGGTGTCAGCAGGGCTGTGATGTCACCAGTCTTGGCTGCAGTGTGGATTAGGTAGGCGGAGGCTGGGAGGCAGTGGGAAGTGGCCGCAGGGACCGTGTGGTCAGTCAGGGTGGTGCCATTTCCTGTGTGCACACTCGGGGGCTCCCGCCACACCAGGGTGCAGACCACCCCGCCGCTCTCCCGCCAGGCAGCCTGCAGGGTGAAGGACTGACATGGCGGGAGGAGTGGCCTTCTCTGTTGGGGCGCCTGTGGCCCCCCAGCTCTTTGATAGGCAGAGCTGCAGGATCTGCCTCCTTGCCCCTGTGAGCCCAGCAGTGATTCTTGTGTCCTCCCCAGATCCATGGCTGCGGTGAGAGGGCTCGCGCCCATGCCCGGCAGAGAATCAGCCGGGAAGGGGTCCTCTATGCCGGCAGCGGGACCAAGAACGGATCCCTGGACCCAGCCAAGAGGGCCCAGCTGCAGAGGAGGCTGGATAAGAAGCTGAGTGAGCTCAGCAACCAGAGGACCAGCCGGAGGAAGGAGAGGGGGACGTGACCGGCCGCATCCTTGCACGCCCCGCGGAGCTCTCTCCATGGTAGCCCAGGGCGCTGGCAGACCATGCTCCGCCTCCACCAGGGCCACAGAGGAGCGGAGGGGCCTATGGGGGAGGAGCGGAGGGCCCTGTTGGGGAAGGTTGGGTTTTTGGACCCCAGGGATAAGCTTTTCCGATGTCACAATGTGGAGGAAAGCACCTGGGGGACAACAGTGCTCGTGCAGGTGGGGCTTGGGAAATGCACGTCCCTTCCCCTTACTCCCCGCCAAAACCCACATCCCAGCCTCTGGATCCTGGGGAAGGTTCCAGTCCCTGGAGAATACCCAGGGCCTCAAACTTGAAGTCACTCCTCCAATGTCTGGGACTTGCCAGCTCAGCCCGTTAGGATGAGGGTGCTGAGAGGAAACAGGAAACAAGACTGCGAATGGCG
</code></pre>

<p>The <code>blastn</code> output from the query above does not contain either of these two Primer-BLAST hits for 'rev_primer' that comprise a potential amplicon from chromosome 11:</p>

<pre><code>product length = 110
Features flanking this product:
   12759 bp at 5' side: neurexin-2-beta isoform X3
   735 bp at 3' side: RAS guanyl-releasing protein 2 isoform X5

Reverse primer  1         CCCAACAGGGCCCTC  15
Template        64726566  AGT............  64726552

Reverse primer  1         CCCAACAGGGCCCTC  15
Template        64726457  ....G..........  64726471
</code></pre>

<p>Why aren't either of these two hits returned in my <code>blastn</code> query? Thanks in advance for helping me with this.</p>

<p><strong>Update:</strong></p>

<p>I noticed in IGV that off-target sequence above is lower-case in the hg38 reference genome. I've tried passing <code>-soft_masking false</code> to blastn, but got the same result.</p>

<p>I thought maybe the problem was due to my using hg38 from UCSC instead of that from NCBI, so I downloaded hg19 from <a href="ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38" rel="nofollow noreferrer">this page</a> and made the blast database from that file. Same result.</p>

<p>I then unmasked the hg38.fna file by converting all lowercase nt to uppercase, as recommended in the NCBI genome assemblies <a href="ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/README.txt" rel="nofollow noreferrer">README.txt</a>, created a BLAST database from that and ran my code again using the unmasked database. Same result, still am not detecting the off-target hit above.</p>

<blast><primer><makeblastdb>
Hits in Primer-BLAST not found with programmatic blastn query
# Hits in Primer-BLAST not found with programmatic blastn query

(Asked by: username_0 on 2019-03-15)

Summary:
I'm writing an automated tool for designing Sanger sequencing primers (other tools don't do exactly what I want) that will be incorporated in a production diagnostic laboratory for on-demand NGS variant verification. I'm stuck on a problem wherein off-target BLAST hits found by NCBI's Primer-BLAST aren't showing up in my programmatic BLAST query of the same genome assembly. 
Background:
After each iteration of design with primer3, I'm using blastn -task blastn to search a local hg38 (patch 12) database for off-target hits. I have logic that does round-robin pairwise comparisons (i.e. compares every hit to every other hit) to determine if each hit pair could result in an amplicon (on same chromosome, opposite-stranded, perfect match at 3'-ends, plus-strand hit is upstream of minus-strand hit, are separated by <= 1000 bp). If my logic finds a viable off-target amplicon, it triggers a total re-design with different design parameters. If no off-targets are found, those primers are finalized, so it's important that I find any and all off-targets.
To make sure that code for viable off-target amplicons is working correctly (that it's triggering re-design or reporting a design failure when appropriate), I'm using NCBI's browser-based Primer-BLAST to verify the absence of viable off-targets, which also uses genome build hg38 patch 12.
The problem:
Primer-BLAST is finding hits that I can't find with my code. Even when I look at every hit returned from the blastn query, I can't find hits that Primer-BLAST is finding. I would like to know if there is some parameter that I am missing, or parameters I could use in addition to those I'm already using, to produce the hits I am getting from Primer-BLAST.
I believe my blast database an my query fasta file are built correctly, because I am getting logical output from the query below.
Programmatic BLAST Query Details:
Query file (query.fasta):
>for_primer
GGATCTGCCTCCTTGCC
>rev_primer
CCCAACAGGGCCCTC

BLASTn Query Statement:
blastn -db blastn -db /path/to/hg38.fasta -evalue 100 -num_threads 8 -outfmt 6 std gaps nident -dust no -word_size 7 -query /path/to/query.fasta

The BLAST database used hg38.fa.gz (which I believe to be patch 12) from here. That was unzipped and processed with makeblastdb.
Primer-BLAST Results:
The Primer-BLAST results are here. To re-create, use the template DNA below, search database = "Genomes for selected organisms...", organism = Homo sapiens, max target size = 1000. The template DNA is:
 GAGGAGTCAGGCCTGATCTGAAGGCAGTGGGGTACCCTCTAGGGTGTCAGCAGGGCTGTGATGTCACCAGTCTTGGCTGCAGTGTGGATTAGGTAGGCGGAGGCTGGGAGGCAGTGGGAAGTGGCCGCAGGGACCGTGTGGTCAGTCAGGGTGGTGCCATTTCCTGTGTGCACACTCGGGGGCTCCCGCCACACCAGGGTGCAGACCACCCCGCCGCTCTCCCGCCAGGCAGCCTGCAGGGTGAAGGACTGACATGGCGGGAGGAGTGGCCTTCTCTGTTGGGGCGCCTGTGGCCCCCCAGCTCTTTGATAGGCAGAGCTGCAGGATCTGCCTCCTTGCCCCTGTGAGCCCAGCAGTGATTCTTGTGTCCTCCCCAGATCCATGGCTGCGGTGAGAGGGCTCGCGCCCATGCCCGGCAGAGAATCAGCCGGGAAGGGGTCCTCTATGCCGGCAGCGGGACCAAGAACGGATCCCTGGACCCAGCCAAGAGGGCCCAGCTGCAGAGGAGGCTGGATAAGAAGCTGAGTGAGCTCAGCAACCAGAGGACCAGCCGGAGGAAGGAGAGGGGGACGTGACCGGCCGCATCCTTGCACGCCCCGCGGAGCTCTCTCCATGGTAGCCCAGGGCGCTGGCAGACCATGCTCCGCCTCCACCAGGGCCACAGAGGAGCGGAGGGGCCTATGGGGGAGGAGCGGAGGGCCCTGTTGGGGAAGGTTGGGTTTTTGGACCCCAGGGATAAGCTTTTCCGATGTCACAATGTGGAGGAAAGCACCTGGGGGACAACAGTGCTCGTGCAGGTGGGGCTTGGGAAATGCACGTCCCTTCCCCTTACTCCCCGCCAAAACCCACATCCCAGCCTCTGGATCCTGGGGAAGGTTCCAGTCCCTGGAGAATACCCAGGGCCTCAAACTTGAAGTCACTCCTCCAATGTCTGGGACTTGCCAGCTCAGCCCGTTAGGATGAGGGTGCTGAGAGGAAACAGGAAACAAGACTGCGAATGGCG

The blastn output from the query above does not contain either of these two Primer-BLAST hits for 'rev_primer' that comprise a potential amplicon from chromosome 11:
product length = 110
Features flanking this product:
   12759 bp at 5' side: neurexin-2-beta isoform X3
   735 bp at 3' side: RAS guanyl-releasing protein 2 isoform X5

Reverse primer  1         CCCAACAGGGCCCTC  15
Template        64726566  AGT............  64726552

Reverse primer  1         CCCAACAGGGCCCTC  15
Template        64726457  ....G..........  64726471

Why aren't either of these two hits returned in my blastn query? Thanks in advance for helping me with this.
Update:
I noticed in IGV that off-target sequence above is lower-case in the hg38 reference genome. I've tried passing -soft_masking false to blastn, but got the same result.
I thought maybe the problem was due to my using hg38 from UCSC instead of that from NCBI, so I downloaded hg19 from this page and made the blast database from that file. Same result.
I then unmasked the hg38.fna file by converting all lowercase nt to uppercase, as recommended in the NCBI genome assemblies README.txt, created a BLAST database from that and ran my code again using the unmasked database. Same result, still am not detecting the off-target hit above.

**Tags:** blast, primer, makeblastdb

---

## 1 answers

### Answer by username_0 on 2019-03-28. Score: 1 | Accepted answer

I've found the problem, and it's way less interesting than I initially thought. Essentially, my "eValue" was way too low. Instead of `-evalue 100`, I'm now using `-evalue 5000` (Primer-BLAST actually uses 30,000).

For some primers, `-evalue 5000` gives ~20,000 results. It's beyond the scope of this question, but I'm now filtering on the query/subject Tm for each hit, perfect match at the query 3'-end, and number of query/subject mismatches. My pairwise
hit comparisons are now detecting plausible off-target amplicons that Primer-BLAST is detecting.

Essentially, I was relying on NCBI's hit ranking system to determine what hits were relevant, but I needed cast a wider net with blastn, and then implement a rational ranking scheme of my own.


---

