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.
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.
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
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:
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.
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.