Hits in Primer-BLAST not found with programmatic blastn query

Question

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.

Answer 1

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.