#!/usr/bin/env python
"""Unit tests for the GenBank database parsers.
"""
from cogent.parse.genbank import parse_locus, parse_single_line, \
indent_splitter, parse_sequence, block_consolidator, parse_organism, \
parse_feature, location_line_tokenizer, parse_simple_location_segment, \
parse_location_line, parse_reference, parse_source, \
Location, LocationList, RichGenbankParser
from cogent.util.unit_test import TestCase, main
__author__ = "Rob Knight"
__copyright__ = "Copyright 2007-2012, The Cogent Project"
__credits__ = ["Rob Knight", "Gavin Huttley"]
__license__ = "GPL"
__version__ = "1.5.3-dev"
__maintainer__ = "Rob Knight"
__email__ = "rob@spot.colorado.edu"
__status__ = "Production"
class GenBankTests(TestCase):
"""Tests of the GenBank main functions."""
def test_parse_locus(self):
"""parse_locus should give correct results on specimen locus lines"""
line = 'LOCUS AF108830 5313 bp mRNA linear PRI 19-MAY-1999'
result = parse_locus(line)
self.assertEqual(len(result), 6)
self.assertEqual(result['locus'], 'AF108830')
self.assertEqual(result['length'], 5313) #note: int, not str
self.assertEqual(result['mol_type'], 'mRNA')
self.assertEqual(result['topology'], 'linear')
self.assertEqual(result['db'], 'PRI')
self.assertEqual(result['date'], '19-MAY-1999')
#should work if some of the fields are missing
line = 'LOCUS AF108830 5313'
result = parse_locus(line)
self.assertEqual(len(result), 2)
self.assertEqual(result['locus'], 'AF108830')
self.assertEqual(result['length'], 5313) #note: int, not str
def test_parse_single_line(self):
"""parse_single_line should split off the label and return the rest"""
line_1 = 'VERSION AF108830.1 GI:4868112\n'
self.assertEqual(parse_single_line(line_1), 'AF108830.1 GI:4868112')
#should work if leading spaces
line_2 = ' VERSION AF108830.1 GI:4868112\n'
self.assertEqual(parse_single_line(line_2), 'AF108830.1 GI:4868112')
def test_indent_splitter(self):
"""indent_splitter should split lines at correct locations"""
#if lines have same indent, should not group together
lines = [
'abc xxx',
'def yyy'
]
self.assertEqual(list(indent_splitter(lines)),\
[[lines[0]], [lines[1]]])
#if second line is indented, should group with first
lines = [
'abc xxx',
' def yyy'
]
self.assertEqual(list(indent_splitter(lines)),\
[[lines[0], lines[1]]])
#if both lines indented but second is more, should group with first
lines = [
' abc xxx',
' def yyy'
]
self.assertEqual(list(indent_splitter(lines)),\
[[lines[0], lines[1]]])
#if both lines indented equally, should not group
lines = [
' abc xxx',
' def yyy'
]
self.assertEqual(list(indent_splitter(lines)), \
[[lines[0]], [lines[1]]])
#for more complex situation, should produce correct grouping
lines = [
' xyz', #0 -
' xxx', #1 -
' yyy', #2
' uuu', #3
' iii', #4
' qaz', #5 -
' wsx', #6 -
' az', #7
' sx', #8
' gb',#9
' bg', #10
' aaa', #11 -
]
self.assertEqual(list(indent_splitter(lines)), \
[[lines[0]], lines[1:5], [lines[5]], lines[6:11], [lines[11]]])
#real example from genbank file
lines = \
"""LOCUS NT_016354 92123751 bp DNA linear CON 29-AUG-2006
DEFINITION Homo sapiens chromosome 4 genomic contig, reference assembly.
ACCESSION NT_016354 NT_006109 NT_006204 NT_006245 NT_006302 NT_006371
NT_006397 NT_016393 NT_016589 NT_016599 NT_016606 NT_022752
NT_022753 NT_022755 NT_022760 NT_022774 NT_022797 NT_022803
NT_022846 NT_022960 NT_025694 NT_028147 NT_029273 NT_030643
NT_030646 NT_030662 NT_031780 NT_031781 NT_031791 NT_034703
NT_034705 NT_037628 NT_037629 NT_079512
VERSION NT_016354.18 GI:88977422
KEYWORDS .
SOURCE Homo sapiens (human)
ORGANISM Homo sapiens
Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
Mammalia; Eutheria; Euarchontoglires; Primates; Haplorrhini;
Catarrhini; Hominidae; Homo.
?
REFERENCE 2 (bases 1 to 92123751)
AUTHORS International Human Genome Sequencing Consortium.
TITLE Finishing the euchromatic sequence of the human genome""".split('\n')
self.assertEqual(list(indent_splitter(lines)), \
[[lines[0]],[lines[1]],lines[2:8],[lines[8]],[lines[9]],lines[10:15],\
[lines[15]], lines[16:]])
def test_parse_sequence(self):
"""parse_sequence should strip bad chars out of sequence lines"""
lines = """
ORIGIN
1 gggagcgcgg cgcgggagcc cgaggctgag actcaccgga ggaagcggcg cgagcgcccc
61 gccatcgtcc \t\t cggctgaagt 123 \ngcagtg \n
121 cctgggctta agcagtcttc45ccacctcagc
//\n\n\n""".split('\n')
result = parse_sequence(lines)
self.assertEqual(result, 'gggagcgcggcgcgggagcccgaggctgagactcaccggaggaagcggcgcgagcgccccgccatcgtcccggctgaagtgcagtgcctgggcttaagcagtcttcccacctcagc')
def test_block_consolidator(self):
"""block_consolidator should join the block together."""
lines = """ ORGANISM Homo sapiens
Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
Mammalia; Eutheria; Euarchontoglires; Primates; Catarrhini;
Hominidae; Homo.""".split('\n')
label, data = block_consolidator(lines)
self.assertEqual(label, 'ORGANISM')
self.assertEqual(data, ['Homo sapiens',
' Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;',
' Mammalia; Eutheria; Euarchontoglires; Primates; Catarrhini;',
' Hominidae; Homo.'])
lines = r"""COMMENT
Contact: Spindel ER
Division of Neuroscience""".splitlines()
label, data = block_consolidator(lines)
self.assertEqual(label, "COMMENT")
self.assertEqual(data, ['', ' Contact: Spindel ER',
' Division of Neuroscience'])
def test_parse_organism(self):
"""parse_organism should return species, taxonomy (up to genus)"""
#note: lines modified to include the following:
# - multiword names
# - multiword names split over a line break
# - periods and other punctuation in names
lines = """ ORGANISM Homo sapiens
Eukaryota; Metazoa; Chordata Craniata; Vertebrata; Euteleostomi;
Mammalia; Eutheria; Euarchontoglires; Primates \t abc. 2.; Catarrhini
Hominidae; Homo.""".split('\n')
species, taxonomy = parse_organism(lines)
self.assertEqual(species, 'Homo sapiens')
self.assertEqual(taxonomy, ['Eukaryota', 'Metazoa', \
'Chordata Craniata', 'Vertebrata', 'Euteleostomi', 'Mammalia', \
'Eutheria', 'Euarchontoglires', 'Primates abc. 2.', \
'Catarrhini Hominidae', 'Homo'])
def test_parse_feature(self):
"""parse_feature should return dict containing annotations of feature"""
example_feature=\
""" CDS complement(join(102262..102647,105026..105217,
106638..106719,152424..152682,243209..243267))
/gene="nad1"
/note="Protein sequence is in conflict with the conceptual
translation; author given translation (not conceptual
translation)
start codon is created by C to U RNA editing"
/codon_start=1
/exception="RNA editing"
/product="NADH dehydrogenase subunit 1"
/protein_id="NP_064011.1"
/db_xref="GI:9838451"
/db_xref="IPI:12345"
/translation="MYIAVPAEILGIILPLLLGVAFLVLAERKVMAFVQRRKGPDVVG
SFGLLQPLADGSKLILKEPISPSSANFSLFRMAPVTTFMLSLVARAVVPFDYGMVLSD
PNIGLLYLFAISSLGVYGIIIAGWSSNSKYAFLGALRSAAQMVPYEVSIGLILITVLI
CVGPRNSSEIVMAQKQIWSGIPLFPVLVMFFISCLAETNRAPFDLPEAERELVAGYNV
EYSSMGSALFFLGEYANMILMSGLCTSLSPGGWPPILDLPISKRIPGSIWFSIKVILF
LFLYIWVRAAFPRYRYDQLMGLGRKVFLPLSLARVVAVSGVLVTFQWLP"""
result = parse_feature(example_feature.split('\n'))
self.assertEqual(result['type'], 'CDS')
self.assertEqual(result['raw_location'], \
['complement(join(102262..102647,105026..105217,', \
' 106638..106719,152424..152682,243209..243267))'])
self.assertEqual(result['gene'], ['nad1'])
self.assertEqual(result['note'], ['Protein sequence is in conflict with the conceptual translation; author given translation (not conceptual translation) start codon is created by C to U RNA editing'])
self.assertEqual(result['codon_start'], ['1'])
self.assertEqual(result['exception'], ['RNA editing'])
self.assertEqual(result['product'], ['NADH dehydrogenase subunit 1'])
self.assertEqual(result['protein_id'],['NP_064011.1'])
self.assertEqual(result['db_xref'], ['GI:9838451','IPI:12345'])
self.assertEqual(result['translation'],['MYIAVPAEILGIILPLLLGVAFLVLAERKVMAFVQRRKGPDVVGSFGLLQPLADGSKLILKEPISPSSANFSLFRMAPVTTFMLSLVARAVVPFDYGMVLSDPNIGLLYLFAISSLGVYGIIIAGWSSNSKYAFLGALRSAAQMVPYEVSIGLILITVLICVGPRNSSEIVMAQKQIWSGIPLFPVLVMFFISCLAETNRAPFDLPEAERELVAGYNVEYSSMGSALFFLGEYANMILMSGLCTSLSPGGWPPILDLPISKRIPGSIWFSIKVILFLFLYIWVRAAFPRYRYDQLMGLGRKVFLPLSLARVVAVSGVLVTFQWLP'])
self.assertEqual(len(result), 11)
short_feature = ['D-loop 15418..16866']
result = parse_feature(short_feature)
self.assertEqual(result['type'], 'D-loop')
self.assertEqual(result['raw_location'], ['15418..16866'])
#can get more than one = in a line
#from AF260826
bad_feature = \
""" tRNA 1173..1238
/note="codon recognized: AUC; Cove score = 16.56"
/product="tRNA-Ile"
/anticodon=(pos:1203..1205,aa:Ile)"""
result = parse_feature(bad_feature.split('\n'))
self.assertEqual(result['note'], \
['codon recognized: AUC; Cove score = 16.56'])
#need not always have an = in a line
#from NC_001807
bad_feature = \
''' mRNA 556
/partial
/citation=[6]
/product="H-strand"'''
result = parse_feature(bad_feature.split('\n'))
self.assertEqual(result['partial'], [''])
def test_location_line_tokenizer(self):
"""location_line_tokenizer should tokenize location lines"""
llt =location_line_tokenizer
self.assertEqual(list(llt(['123..456'])), ['123..456'])
self.assertEqual(list(llt(['complement(123..456)'])), \
['complement(', '123..456', ')'])
self.assertEqual(list(llt(['join(1..2,3..4)'])), \
['join(', '1..2', ',', '3..4', ')'])
self.assertEqual(list(llt([\
'join(complement(1..2, join(complement( 3..4),',\
'\n5..6), 7..8\t))'])),\
['join(','complement(','1..2',',','join(','complement(','3..4',\
')', ',', '5..6',')',',','7..8',')',')'])
def test_parse_simple_location_segment(self):
"""parse_simple_location_segment should parse simple segments"""
lsp = parse_simple_location_segment
l = lsp('37')
self.assertEqual(l._data, 37)
self.assertEqual(str(l), '37')
self.assertEqual(l.Strand, 1)
l = lsp('40..50')
first, second = l._data
self.assertEqual(first._data, 40)
self.assertEqual(second._data, 50)
self.assertEqual(str(l), '40..50')
self.assertEqual(l.Strand, 1)
#should handle ambiguous starts and ends
l = lsp('>37')
self.assertEqual(l._data, 37)
self.assertEqual(str(l), '>37')
l = lsp('<37')
self.assertEqual(l._data, 37)
self.assertEqual(str(l), '<37')
l = lsp('<37..>42')
first, second = l._data
self.assertEqual(first._data, 37)
self.assertEqual(second._data, 42)
self.assertEqual(str(first), '<37')
self.assertEqual(str(second), '>42')
self.assertEqual(str(l), '<37..>42')
def test_parse_location_line(self):
"""parse_location_line should give correct list of location objects"""
llt = location_line_tokenizer
r = parse_location_line(llt(['123..456']))
self.assertEqual(str(r), '123..456')
r = parse_location_line(llt(['complement(123..456)']))
self.assertEqual(str(r), 'complement(123..456)')
r = parse_location_line(llt(['complement(123..456, 345..678)']))
self.assertEqual(str(r), \
'join(complement(345..678),complement(123..456))')
r = parse_location_line(llt(['complement(join(123..456, 345..678))']))
self.assertEqual(str(r), \
'join(complement(345..678),complement(123..456))')
r = parse_location_line(\
llt(['join(complement(123..456), complement(345..678))']))
self.assertEqual(str(r), \
'join(complement(123..456),complement(345..678))')
#try some nested joins and complements
r = parse_location_line(llt(\
['complement(join(1..2,3..4,complement(5..6),',
'join(7..8,complement(9..10))))']))
self.assertEqual(str(r), \
'join(9..10,complement(7..8),5..6,complement(3..4),complement(1..2))')
def test_parse_reference(self):
"""parse_reference should give correct fields"""
r = \
"""REFERENCE 2 (bases 1 to 2587)
AUTHORS Janzen,D.M. and Geballe,A.P.
TITLE The effect of eukaryotic release factor depletion on translation
termination in human cell lines
JOURNAL (er) Nucleic Acids Res. 32 (15), 4491-4502 (2004)
PUBMED 15326224"""
result = parse_reference(r.split('\n'))
self.assertEqual(len(result), 5)
self.assertEqual(result['reference'], '2 (bases 1 to 2587)')
self.assertEqual(result['authors'], 'Janzen,D.M. and Geballe,A.P.')
self.assertEqual(result['title'], \
'The effect of eukaryotic release factor depletion ' + \
'on translation termination in human cell lines')
self.assertEqual(result['journal'], \
'(er) Nucleic Acids Res. 32 (15), 4491-4502 (2004)')
self.assertEqual(result['pubmed'], '15326224')
def test_parse_source(self):
"""parse_source should split into source and organism"""
s = \
"""SOURCE African elephant.
ORGANISM Mitochondrion Loxodonta africana
Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
Mammalia; Eutheria; Proboscidea; Elephantidae; Loxodonta.""".split('\n')
r = parse_source(s)
self.assertEqual(len(r), 3)
self.assertEqual(r['source'], 'African elephant.')
self.assertEqual(r['species'], 'Mitochondrion Loxodonta africana')
self.assertEqual(r['taxonomy'], ['Eukaryota','Metazoa', 'Chordata',\
'Craniata', 'Vertebrata', 'Euteleostomi', 'Mammalia',\
'Eutheria', 'Proboscidea', 'Elephantidae', 'Loxodonta'])
def test_rich_parser(self):
"""correctly constructs +/- strand features"""
# a custom annotation function
from cogent.core.annotation import Feature
def add_annotation(seq, feature, spans):
if feature['type'] != "CDS": return
name = feature['locus_tag'][0]
seq.addAnnotation(Feature, "CDS", name, spans)
parser = RichGenbankParser(open('data/annotated_seq.gb'),
add_annotation=add_annotation)
seq = [s for l, s in parser][0]
cds = dict([(f.Name, f) for f in seq.getAnnotationsMatching('CDS')])
expects = {
'CNA00110': 'MAGYDARYGNPLDPMSGGRPSPPETSQQDAYEYSKHGSSSGYLGQLPLGAD'\
'SAQAETASALRTLFGEGADVQALQEPPNQINTLAEGAAVAETGGVLGGDTTRSDNEALAIDPSL'\
'SEQAAPAPKDSTETPDDRSRSPSSGNHHHHHPAVKRKATSRAGMLARGGACEFCKRRKLKCSAEL'\
'PACANCVKSGKECVYAQKKQRSRVKVLEDRLQELEKRLEQGQAGAASASGGDSGAHAASSVYTAP'\
'SLGSGGGSELTVEQTLVHNVDPSLLPPSEYDEAFILHDFDSFADMRKQETQLEPDLMTLADAAAA'\
'DTPAAAAAETNDPWAKMSPEEIVKEIIKVATGGKGEGERIISHLVQTYMNSTVNTWHPLVIPPMD'\
'LVSRVSRTTPDPIHPTLLLSLIPALLPLSPIQSLRHPAIPLLLLPHARAHSVQAITQSDPRVLDT'\
'IIAGVSRAYSFFNEAKNIDGWVDCVAATSLVRAAGLTKQGGVGERFVPEDRVPAERLAKRRREAG'\
'LRALMHKGAIVPPPESWYQFGQRVNLFWTSYICDRAAAIGWGWPSSYNDEDITTPWPKDDYKSVQ'\
'ALLDDTTIHTFLSPLAPAPAPATPDSDLCAQAKSITLLYHAQRLLDSPPELSTPEKTHRLLGLTE'\
'GYMESLEKMRGPRMRAGKLSSVWMILYTTIAVLHSKDGFDKCDPDGADQVSITRVVAAADKVLEL'\
'VSAVQNTGDTHLSSCDVISSVLFLHLARLMIQYTNRLRLRVQDSALVSTLRAKTESFKRALIDQG'\
'ERLVFAQVAAQMLENYHVGAEWKAGEWERADGGDWRGV',
'CNA00120': 'MDFSQFNGAEQAHMSKVIEKKQMQDFMRLYSGLVEKCFNACAQD'\
'FTSKALTTNETTCVQNCTDKFLKHSERVGARFAEHNAGMLSPYGAASLMASQSKCRAP'\
'DSNGLGVFCKWRRIKSTVVLYNHLACIKQMDNRF'}
for locus in cds:
got = cds[locus].getSlice().\
withoutTerminalStopCodon().getTranslation()
self.assertEqual(str(got), expects[locus])
class LocationTests(TestCase):
"""Tests of the Location class."""
def test_init(self):
"""Location should init with 1 or 2 values, plus params."""
l = Location(37)
self.assertEqual(str(l), '37')
l = Location(37, Ambiguity = '>')
self.assertEqual(str(l), '>37')
l = Location(37, Ambiguity='<')
self.assertEqual(str(l), '<37')
l = Location(37, Accession='AB123')
self.assertEqual(str(l), 'AB123:37')
l = Location(37, Accession='AB123', Db='Kegg')
self.assertEqual(str(l), 'Kegg::AB123:37')
l1 = Location(37)
l2 = Location(42)
l = Location([l1,l2])
self.assertEqual(str(l), '37..42')
l3 = Location([l1,l2], IsBounds=True)
self.assertEqual(str(l3), '(37.42)')
l4 = Location([l1,l2], IsBetween=True)
self.assertEqual(str(l4), '37^42')
l5 = Location([l4,l3])
self.assertEqual(str(l5), '37^42..(37.42)')
l5 = Location([l4,l3], Strand=-1)
self.assertEqual(str(l5), 'complement(37^42..(37.42))')
class LocationListTests(TestCase):
"""Tests of the LocationList class."""
def test_extract(self):
"""LocationList extract should return correct sequence"""
l = Location(3)
l2_a = Location(5)
l2_b = Location(7)
l2 = Location([l2_a,l2_b], Strand=-1)
l3_a = Location(10)
l3_b = Location(12)
l3 = Location([l3_a, l3_b])
ll = LocationList([l, l2, l3])
s = ll.extract('ACGTGCAGTCAGTAGCAT')
# 123456789012345678
self.assertEqual(s, 'G'+'TGC'+'CAG')
#check a case where it wraps around
l5_a = Location(16)
l5_b = Location(4)
l5 = Location([l5_a,l5_b])
ll = LocationList([l5])
s = ll.extract('ACGTGCAGTCAGTAGCAT')
self.assertEqual(s, 'CATACGT')
if __name__ == '__main__':
from sys import argv
if len(argv) > 2 and argv[1] == 'x':
filename = argv[2]
lines = open(filename)
for i in indent_splitter(lines):
print '******'
print i[0]
for j in indent_splitter(i[1:]):
print '?????'
for line in j:
print line
else:
main()
