Extraction of LRRK2 gene coding regions

LRRK2 gene includes 51 coding exons (7784 bp including two intron nucleotides near exon-intron junctions). To design NGS-panel for sequencing it, first, you need to extract genome coordinates of regions studied. To do it, create text file with one word in it - LRRK2 and run script getGeneRegions.py:

python3 getGeneRegions.py \
-glf LRRK2.txt \
-ref <path to the directory with reference genome GenBank-files> \
-wgref <path to the reference genome FASTA-file> \
-rf LRRK2.regions.csv

Open the output file in Microsoft Excel (or other similar programs) and change number of multiplexes for each region in 5th column from 1 to necessary number of multiplex reactions. This number can be estimated from the following facts: (1) total number of primer pairs (about 50-70 primer pairs can be joined together); (2) how many overlapping regions (overlapped primer pairs can not be joined); (3) size of amplicons (larger amplicons - less overlappings). Initially, we tried to sort primer pairs to 3 multiplex reactions, so change 1 in the 5th column from 1 to "1,2,3" (without quotes) and save it.

Running primer design

Default settings are for DNA samples from formalin-fixed paraffin-embedded (FFPE) tissue specimens. To study germline mutations, like in this case, we can use amplicons of about 150 bp. So, we will use Minimal amplicon length as 130 bp, optimal and maximal - as 150 bp. Also, you can decrease number of primer pairs designed for 1 locus (-primernum1) to 5 (then the process will be faster) and set number of nucleotides from 3'-end to check for overlapping with SNP (-nucs) to 10. We can leave other values by default. So, the command will be the following:

python3 NGS_primerplex.py \
-regions LRRK2.regions.csv \
-ref <path to reference genome FASTA-file> \
-ad1 cgaccacgctcttgaagca \
-ad2 gggttccctaagggttgga \
-minampllen 130 \
-optampllen 150 \
-maxampllen 150 \
-blast \
-snps \
-nucs 10 \
-dbsnp <path to VCF-file with SNPs> \
-th <number of threads that you computer allows using> \
-primernum1 5

Initial primer design will take about 1-3 hours. When it finished, you likely got that for 241 regions no primers could be designed. So, we need to use less strict parameters, e.g. we can increase maximal length of primers. Now it is 28 nucleotides, try to increase it to 35, and start design again with choosing draft primers (primers that program could design):

python3 NGS_primerplex.py \
-regions LRRK2.regions.csv \
-ref <path to reference genome FASTA-file> \
-ad1 cgaccacgctcttgaagca \
-ad2 gggttccctaagggttgga \
-minampllen 130 \
-optampllen 150 \
-maxampllen 150 \
-blast \
-snps \
-nucs 10 \
-dbsnp <path to VCF-file with SNPs> \
-th <number of threads that you computer allows using> \
-primernum1 5 \
-maxprimerlen 35 \
-draft LRRK2.regions_NGS_primerplex_all_draft_primers.xls

After 10 minutes, you likely got that for 24 regions primers couldn't be designed. Try to increase maximal primer melt from 68 to 72 and minimal primer GC-content from 20 to 15. Note, that you don't have to change draft-file because it is updated automatically:

python3 NGS_primerplex.py \
-regions LRRK2.regions.csv \
-ref <path to reference genome FASTA-file> \
-ad1 cgaccacgctcttgaagca \
-ad2 gggttccctaagggttgga \
-minampllen 130 \
-optampllen 150 \
-maxampllen 150 \
-blast \
-snps \
-nucs 10 \
-dbsnp <path to VCF-file with SNPs> \
-th <number of threads that you computer allows using> \
-primernum1 5 \
-maxprimerlen 35 \
-maxprimermelt 70 \
-minprimergc 15 \
-draft LRRK2.regions_NGS_primerplex_all_draft_primers.xls

After 1 minute you probably got that for 18 regions primers couldn't be designed. Then we can try to make less strict parameters for secondary structures for one primer pair (for dimers without hybridized 3'-end from 35 to 40, for hairpins from 40 to 45):

python3 NGS_primerplex.py \
-regions LRRK2.regions.csv \
-ref <path to reference genome FASTA-file> \
-ad1 cgaccacgctcttgaagca \
-ad2 gggttccctaagggttgga \
-minampllen 130 \
-optampllen 150 \
-maxampllen 150 \
-blast \
-snps \
-nucs 10 \
-dbsnp <path to VCF-file with SNPs> \
-th <number of threads that you computer allows using> \
-primernum1 5 \
-maxprimerlen 35 \
-maxprimermelt 70 \
-minprimergc 15 \
-draft LRRK2.regions_NGS_primerplex_all_draft_primers.xls \
-maxprimercomplanyth 40 \
-maxprimerhairpinth 45

After 1 minute more we got only 3 regions without primers. We can try to decrease minimal primer melting temperature from 60 to 59:

python3 NGS_primerplex.py \
-regions LRRK2.regions.csv \
-ref <path to reference genome FASTA-file> \
-ad1 cgaccacgctcttgaagca \
-ad2 gggttccctaagggttgga \
-minampllen 130 \
-optampllen 150 \
-maxampllen 150 \
-blast \
-snps \
-nucs 10 \
-dbsnp <path to VCF-file with SNPs> \
-th <number of threads that you computer allows using> \
-primernum1 5 \
-maxprimerlen 35 \
-maxprimermelt 70 \
-minprimergc 15 \
-draft LRRK2.regions_NGS_primerplex_all_draft_primers.xls \
-maxprimercomplanyth 40 \
-maxprimerhairpinth 45 \
-minprimermelt 59

Now program goes through this step. But after checking for covering SNPs we get some regions for which all primers were removed due to overlapping with SNP. So, we need to change draft-file to file after SNPs and increase number of primer pairs designed for 1 locus (-primernum1). This will make NGS-PrimerPlex to design more primers from which maybe we will be able to separate primers not covering SNPs. Now we can return our strict parameters to default:

python3 NGS_primerplex.py \
-regions LRRK2.regions.csv \
-ref <path to reference genome FASTA-file> \
-ad1 cgaccacgctcttgaagca \
-ad2 gggttccctaagggttgga \
-minampllen 130 \
-optampllen 150 \
-maxampllen 150 \
-blast \
-snps \
-nucs 10 \
-dbsnp <path to VCF-file with SNPs> \
-th <number of threads that you computer allows using> \
-primernum1 100 \
-draft LRRK2.regions_NGS_primerplex_all_draft_primers_after_SNPs.xls

For 29 regions primers couldn't be designed. Returning less strict parameters...

python3 NGS_primerplex.py \
-regions LRRK2.regions.csv \
-ref <path to reference genome FASTA-file> \
-ad1 cgaccacgctcttgaagca \
-ad2 gggttccctaagggttgga \
-minampllen 130 \
-optampllen 150 \
-maxampllen 150 \
-blast \
-snps \
-nucs 10 \
-dbsnp <path to VCF-file with SNPs> \
-th <number of threads that you computer allows using> \
-primernum1 100 \
-draft LRRK2.regions_NGS_primerplex_all_draft_primers_after_SNPs.xls
-maxprimerlen 35 \
-maxprimermelt 70 \
-minprimergc 15

Now, for 36 regions no primers left after filtering by SNPs (versus 154 previously). So, increase -primernum1 (e.g. to 200) and start again:

python3 NGS_primerplex.py \
-regions LRRK2.regions.csv \
-ref <path to reference genome FASTA-file> \
-ad1 cgaccacgctcttgaagca \
-ad2 gggttccctaagggttgga \
-minampllen 130 \
-optampllen 150 \
-maxampllen 150 \
-blast \
-snps \
-nucs 10 \
-dbsnp <path to VCF-file with SNPs> \
-th <number of threads that you computer allows using> \
-primernum1 200 \
-draft LRRK2.regions_NGS_primerplex_all_draft_primers_after_SNPs.xls
-maxprimerlen 35 \
-maxprimermelt 70 \
-minprimergc 15

By such subsequent changing number of primer pairs for 1 locus (e.g. 300 and then 400) and choosing less strict parameters (if it is necessary), we got this step completed. But during combining different primer pairs to combinations we got new error: "Too low value of maximal overlap (-maxoverlap) or of initially designed primers (-primernum): 50 and 200. Try to increase one of them". This means that some neighboring amplicons overlap with each other more than value defined (more than 50). We can try to increase it to e.g. 100 and run again.

After that, we got that some primer pairs (13-16) couldn't be sorted to any of multiplex reactions. In this case, when so many primer pairs are unsorted, we can increase number of primer pools or use less strict parameters while checking primers for interactions and non-target product formation (-minmultdimerdg1 and -minmultdimerdg2). First, we will try to decrease minimal acceptable value of free energy of primer dimer formation without hybridized 3'-end (-minmultdimerdg2) from -10 to -12:

python3 NGS_primerplex.py \
-regions LRRK2.regions.csv \
-ref <path to reference genome FASTA-file> \
-ad1 cgaccacgctcttgaagca \
-ad2 gggttccctaagggttgga \
-minampllen 130 \
-optampllen 150 \
-maxampllen 150 \
-blast \
-snps \
-nucs 10 \
-dbsnp <path to VCF-file with SNPs> \
-th <number of threads that you computer allows using> \
-primernum1 200 \
-draft LRRK2.regions_NGS_primerplex_all_draft_primers_after_SNPs.xls
-maxprimerlen 35 \
-maxprimermelt 70 \
-minprimergc 15 \
-maxoverlap 100 \
-minmultdimerdg2 -12