Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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
#!/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()
#!/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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.