Script for performing some basic testing on Illumina amplicon sequencing primers as described in: http://www.nature.com/ismej/journal/v6/n8/full/ismej20128a.html

check_illumina_barcodes.py

#!/usr/bin/env python
# File created on 01 Dec 2011
from __future__ import division

__author__ = "Greg Caporaso"
__copyright__ = "Copyright 2011, The QIIME project"
__credits__ = ["Greg Caporaso"]
__license__ = "GPL"
__version__ = "1.3.0-dev"
__maintainer__ = "Greg Caporaso"
__email__ = "gregcaporaso@gmail.com"
__status__ = "Development"
 

from cogent import DNA
from qiime.util import parse_command_line_parameters, make_option

script_info = {}
script_info['brief_description'] = ""
script_info['script_description'] = ""
script_info['script_usage'] = [("","","""%prog --pcr_f "AATGATACGGCGACCACCGAGATCTACAC TATGGTAATT GT GTGCCAGCMGCCGCGGTAA" --pcr_r "CAAGCAGAAGACGGCATACGAGAT TCCCTTGTCTCC AGTCAGTCAG CC GGACTACHVGGGTWTCTAAT" --seq_r1 "TATGGTAATT GT GTGCCAGCMGCCGCGGTAA" --seq_r2 "AGTCAGTCAG CC GGACTACHVGGGTWTCTAAT" --seq_i "ATTAGAWACCCBDGTAGTCC GG CTGACTGACT" """)]
script_info['output_description']= ""
script_info['required_options'] = [\
 # Example required option
 make_option('--pcr_f',help='Forward primer PCR sequence'),
 make_option('--pcr_r',help='Reverse primer PCR sequence'),
 make_option('--seq_r1',help='Read 1 sequencing primer'),
 make_option('--seq_r2',help='Read 2 sequencing primer'),
 make_option('--seq_i',help='Index sequencing primer'),
]
script_info['optional_options'] = []
script_info['version'] = __version__

def check_read_sequencing_primers(pcr,seq):
    error_messages = []
    pcr = pcr.split()[-3:]
    seq = seq.split()
    if pcr[0] != seq[0]:
        error_messages.append("Non-matching pad regions: %s != %s." % (pcr[0], seq[0]))
    if pcr[1] != seq[1]:
        error_messages.append("Non-matching linker regions: %s != %s." % (pcr[1], seq[1]))
    if pcr[2] != seq[2]:
        error_messages.append("Non-matching primer: %s != %s." % (pcr[2], seq[2]))
    return error_messages

def check_index_sequencing_primers(pcr,seq):
    error_messages = []
    pcr = pcr.split()[-3:]
    seq = DNA.rc(seq).split()
    if pcr[0] != seq[0]:
        error_messages.append("Non-rc pad regions: %s != %s." % (pcr[0], seq[0]))
    if pcr[1] != seq[1]:
        error_messages.append("Non-rc linker regions: %s != %s." % (pcr[1], seq[1]))
    if pcr[2] != seq[2]:
        error_messages.append("Non-rc primer: %s != %s." % (pcr[2], seq[2]))
    return error_messages


def main():
    option_parser, opts, args =\
       parse_command_line_parameters(**script_info)

    print "Read1/Forward PCR primer check:"
    error_messages = check_read_sequencing_primers(opts.pcr_f,opts.seq_r1)
    if error_messages:
        print "There are errors in your primer:"
        print ' ' + '\n '.join(error_messages)
    else:
        print "No errors."
    
    print ''
    print "Read2/Reverse PCR primer check:"
    error_messages = check_read_sequencing_primers(opts.pcr_r,opts.seq_r2)
    if error_messages:
        print "There are errors in your primer:"
        print ' ' + '\n '.join(error_messages)
    else:
        print "No errors."

    print ''
    print "Index read/Reverse PCR primer check:"
    error_messages = check_index_sequencing_primers(opts.pcr_r,opts.seq_i)
    if error_messages:
        print "There are errors in your primer:"
        print ' ' + '\n '.join(error_messages)
    else:
        print "No errors."


if __name__ == "__main__":
    main()

test_check_illumina_primers.py

#!/usr/bin/env python
# File created on 01 Dec 2011
from __future__ import division

__author__ = "Greg Caporaso"
__copyright__ = "Copyright 2011, The QIIME project"
__credits__ = ["Greg Caporaso"]
__license__ = "GPL"
__version__ = "1.3.0-dev"
__maintainer__ = "Greg Caporaso"
__email__ = "gregcaporaso@gmail.com"
__status__ = "Development"
 

from cogent.util.unit_test import TestCase, main
from check_illumina_primers import (check_read_sequencing_primers,
                                    check_index_sequencing_primers)

class CheckIlluminaPrimers(TestCase):
    """ """
    def setUp(self):
        """ """
        self.pcr_f = "AATGATACGGCGACCACCGAGATCTACAC TATGGTAATT GT GTGCCAGCMGCCGCGGTAA"
        self.pcr_r = "CAAGCAGAAGACGGCATACGAGAT TCCCTTGTCTCC AGTCAGTCAG CC GGACTACHVGGGTWTCTAAT"
        self.seq_r1 = "TATGGTAATT GT GTGCCAGCMGCCGCGGTAA"
        self.seq_r2 = "AGTCAGTCAG CC GGACTACHVGGGTWTCTAAT"
        self.seq_i = "ATTAGAWACCCBDGTAGTCC GG CTGACTGACT"
    
    def test_check_read_sequencing_primers_f(self):
        """forward/read1 works as expected"""
        self.assertEqual(check_read_sequencing_primers(self.pcr_f,self.seq_r1),[])
    
    def test_check_read_sequencing_primers_r(self):
        """reverse/read2 works as expected"""
        self.assertEqual(check_read_sequencing_primers(self.pcr_r,self.seq_r2),[])
    
    def test_check_read_sequencing_primers_detects_errors(self):
        """pad error detected in read1/read2 primer check"""
        # error in pad
        seq_r1 = "TATGCTAATT GT GTGCCAGCMGCCGCGGTAA"
        self.assertEqual(len(check_read_sequencing_primers(self.pcr_f,seq_r1)),
                         1)
        # error in linker
        seq_r1 = "TATGGTAATT GA GTGCCAGCMGCCGCGGTAA"
        self.assertEqual(len(check_read_sequencing_primers(self.pcr_f,seq_r1)),
                         1)
        # error in primer
        seq_r1 = "TATGGTAATT GT GTGCCAGCMGCCGGGGTAA"
        self.assertEqual(len(check_read_sequencing_primers(self.pcr_f,seq_r1)),
                         1)
        # error in all three
        seq_r1 = "TATGGCAATT GA GTGCCAGCMGCCGGGGTAA"
        self.assertEqual(len(check_read_sequencing_primers(self.pcr_f,seq_r1)),
                         3)
    
    def test_check_index_sequencing_primers(self):
        """reverse/index functions as expected"""
        self.assertEqual(check_index_sequencing_primers(self.pcr_r,self.seq_i),[])

    def test_check_index_sequencing_primers_detects_errors(self):
        """reverse/index detects errors"""
        # error in pad
        seq_i = "ATTAGAWACCCBDGTAGTCC GG CTGACCGACT"
        self.assertEqual(len(check_index_sequencing_primers(self.pcr_r,seq_i)),
                         1)
        # error in linker 
        seq_i = "ATTAGAWACCCBDGTAGTCC AG CTGACTGACT"
        self.assertEqual(len(check_index_sequencing_primers(self.pcr_r,seq_i)),
                         1)
        # error in primer
        seq_i = "ATTAGAWACCCBDGTAGTCC GG CTGACAGACT"
        self.assertEqual(len(check_index_sequencing_primers(self.pcr_r,seq_i)),
                         1)
        # error in all three
        seq_i = "ATTACAWACCCBDGTAGTCC GC CTGTCTGACT"
        self.assertEqual(len(check_index_sequencing_primers(self.pcr_r,seq_i)),
                         3)


if __name__ == "__main__":
    main()