1
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.