text: user3564: 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.
  
  user3564: 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.
  
original_text: 3564: <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>
   3564: <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>
  