t/SeqIO/embl.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7         use lib '.';
   8         use Bio::Root::Test;
   9
  10         test_begin(-tests => 69);
  11
  12     use_ok('Bio::SeqIO::embl');
  13 }
  14
  15 my $verbose = test_debug();
  16
  17 my $ast = Bio::SeqIO->new( -format => 'embl',
  18                            -verbose => $verbose,
  19                            -file => test_input_file('roa1.dat'));
  20 $ast->verbose($verbose);
  21 my $as = $ast->next_seq();
  22 ok defined $as->seq;
  23 is($as->display_id, 'HSHNCPA1');
  24 is($as->accession_number, 'X79536');
  25 is($as->seq_version, 1);
  26 is($as->version, 1);
  27 is($as->desc, 'H.sapiens mRNA for hnRNPcore protein A1');
  28 is($as->molecule, 'RNA');
  29 is($as->alphabet, 'rna');
  30 is(scalar $as->all_SeqFeatures(), 4);
  31 is($as->length, 1198);
  32 is($as->species->binomial(), 'Homo sapiens');
  33
  34 # EMBL Release 87 changes (8-17-06)
  35
  36 $ast = Bio::SeqIO->new( -format => 'embl',
  37                            -verbose => $verbose,
  38                            -file => test_input_file('roa1_v2.dat'));
  39 $ast->verbose($verbose);
  40 $as = $ast->next_seq();
  41 ok defined $as->seq;
  42 # accession # same as display name now
  43 is($as->display_id, 'X79536');
  44 is($as->accession_number, 'X79536');
  45 is($as->seq_version, 1);
  46 is($as->version, 1);
  47 is($as->desc, 'H.sapiens mRNA for hnRNPcore protein A1');
  48 # mRNA instead of RNA
  49 is($as->molecule, 'mRNA');
  50 is($as->alphabet, 'rna');
  51 is(scalar $as->all_SeqFeatures(), 4);
  52 is($as->length, 1198);
  53 is($as->species->binomial(), 'Homo sapiens');
  54
  55 my $ent = Bio::SeqIO->new( -file => test_input_file('test.embl'),
  56                            -format => 'embl');
  57 my $seq = $ent->next_seq();
  58
  59 is(defined $seq->seq(), 1,
  60    'success reading Embl with ^ location and badly split double quotes');
  61 is(scalar $seq->annotation->get_Annotations('reference'), 3);
  62
  63 my $out_file = test_output_file();
  64 my $out = Bio::SeqIO->new(-file=> ">$out_file",
  65                           -format => 'embl');
  66 is($out->write_seq($seq),1,
  67    'success writing Embl format with ^ < and > locations');
  68
  69 # embl with no FT
  70 $ent = Bio::SeqIO->new( -file => test_input_file('test.embl'),
  71                         -format => 'embl');
  72 $seq = $ent->next_seq();
  73
  74 ok($seq);
  75 is(lc($seq->subseq(1,10)),'gatcagtaga');
  76 is($seq->length, 4870);
  77
  78 # embl with no FH
  79 my $noFH = Bio::SeqIO->new(-file => test_input_file('no_FH.embl'),
  80                         -format => 'embl');
  81 is(scalar($noFH->next_seq->get_SeqFeatures), 4);
  82
  83 # bug 1571
  84 $ent = Bio::SeqIO->new(-format => 'embl',
  85                        -file   => test_input_file('test.embl2sq'));
  86 is($ent->next_seq->length,4877);
  87
  88 # embl repbase
  89 $ent = Bio::SeqIO->new(-file => test_input_file('BEL16-LTR_AG.embl'), -format => 'embl');
  90 $seq = $ent->next_seq;
  91 is($seq->display_id,'BEL16-LTR_AG');
  92
  93 # test secondary accessions in EMBL (bug #1332)
  94 my $seqio = Bio::SeqIO->new(-format => 'embl',
  95                            -file => test_input_file('ECAPAH02.embl'));
  96 $seq = $seqio->next_seq;
  97 is($seq->accession_number, 'D10483');
  98 is($seq->seq_version, 2);
  99 my @accs = $seq->get_secondary_accessions();
 100 is($accs[0], 'J01597');
 101 is($accs[-1], 'X56742');
 102
 103 ### TPA TESTS - Thanks to Richard Adams ###
 104 # test Third Party Annotation entries in EMBL/Gb format
 105 # to ensure compatability with parsers.
 106 my $str = Bio::SeqIO->new(-format =>'embl',
 107                          -file => test_input_file('BN000066-tpa.embl'));
 108 $seq = $str->next_seq;
 109 ok(defined $seq);
 110 is($seq->accession_number, 'BN000066');
 111 is($seq->alphabet, 'dna');
 112 is($seq->display_id, 'AGA000066');
 113 is($seq->length, 5195);
 114 is($seq->division, 'INV');
 115 is($seq->get_dates, 2);
 116 is($seq->keywords, 'acetylcholinesterase; achE1 gene; Third Party Annotation; TPA.');
 117 is($seq->seq_version, 1);
 118 is($seq->feature_count, 15);
 119
 120 my $spec_obj = $seq->species;
 121 is ($spec_obj->common_name, 'African malaria mosquito');
 122 is ($spec_obj->species, 'gambiae');
 123 is ($spec_obj->genus, 'Anopheles');
 124 is ($spec_obj->binomial, 'Anopheles gambiae');
 125
 126 my $ac = $seq->annotation;
 127 my $reference =  ($ac->get_Annotations('reference') )[1];
 128 is ($reference->title,'"A novel acetylcholinesterase gene in mosquitoes codes for the insecticide target and is non-homologous to the ace gene in Drosophila"');
 129 is ($reference->authors,'Weill M., Fort P., Berthomi eu A., Dubois M.P., Pasteur N., Raymond M.');
 130 my $cmmnt =  ($ac->get_Annotations('comment') )[0];
 131 is($cmmnt->text, 'see also AJ488492 for achE-1 from Kisumu strain Third Party Annotation Database: This TPA record uses Anopheles gambiae trace archive data (http://trace.ensembl.org) ');
 132
 133
 134 $ent = Bio::SeqIO->new( -file => test_input_file('test.embl'),
 135                         -format => 'embl');
 136 $ent->verbose($verbose);
 137 $seq = $ent->next_seq();
 138 my $species = $seq->species();
 139 my @cl = $species->classification();
 140 is( $cl[3] ne $species->genus(), 1, 'genus duplication test');
 141 $ent->close();
 142
 143 #
 144 ## read-write - test embl writing of a PrimarySeq
 145 #
 146 my $primaryseq = Bio::PrimarySeq->new( -seq => 'AGAGAGAGATA',
 147                                       -id  => 'myid',
 148                                       -desc => 'mydescr',
 149                                       -alphabet => 'DNA',
 150                                       -accession_number => 'myaccession');
 151
 152
 153
 154 $verbose = -1 unless $ENV{'BIOPERLDEBUG'};  # silence warnings unless we are debuggin
 155
 156 my $embl = Bio::SeqIO->new(-format => 'embl',
 157                           -verbose => $verbose,
 158                           -file => ">$out_file");
 159
 160 ok($embl->write_seq($primaryseq));
 161
 162 # this should generate a warning
 163 my $scalar = "test";
 164 eval {
 165         $embl->write_seq($scalar);
 166 };
 167 ok ($@);
 168
 169 # CDS records
 170 # (which have nonstandard 'PA' and 'OX' tags)
 171 # see http://bioperl.org/pipermail/bioperl-l/2009-February/029252.html
 172 # and the rest of that thread
 173 my $cds_file = Bio::SeqIO->new(-format =>'embl',
 174                                -file => test_input_file('cds_sample.embl'));
 175 my $cds_seq = $cds_file->next_seq;
 176 ok(defined $cds_seq);
 177 is($cds_seq->display_id, 'EAL24309');
 178 is($cds_seq->accession_number, 'CH236947.1', 'CDS - accession on PA line');
 179 is($cds_seq->alphabet, 'dna');
 180 is($cds_seq->length, 192);
 181 is($cds_seq->species->binomial(), 'Homo sapiens');
 182 is($cds_seq->seq_version, 1);
 183 is($cds_seq->feature_count, 2);
 184 my $cds_annot = $cds_seq->annotation;
 185 ok(defined $cds_annot);
 186 my $cds_dblink = ($cds_annot->get_Annotations('dblink'))[0];
 187 ok(defined $cds_dblink);
 188 is($cds_dblink->tagname, 'dblink', 'CDS - OX tagname');
 189 is($cds_dblink->database, 'NCBI_TaxID', 'CDS - OX database');
 190 is($cds_dblink->primary_id, '9606', 'CDS - OX primary_id');
 191