[Bio] / FigKernelPackages / SeedUtils.pm Repository:
ViewVC logotype

Diff of /FigKernelPackages/SeedUtils.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.38, Wed Jun 16 20:21:50 2010 UTC revision 1.39, Wed Jul 21 15:59:36 2010 UTC
# Line 40  Line 40 
40                       rev_comp genome_of min max sims verify_dir between translate                       rev_comp genome_of min max sims verify_dir between translate
41                       standard_genetic_code parse_location roles_of_function                       standard_genetic_code parse_location roles_of_function
42                       strip_ec location_string location_cmp strand_of by_fig_id                       strip_ec location_string location_cmp strand_of by_fig_id
43                       verify_db bbh_data id_url);                       verify_db bbh_data id_url validate_fasta_file);
44    
45  =head1 SEED Utility Methods  =head1 SEED Utility Methods
46    
# Line 1509  Line 1509 
1509      }      }
1510  }  }
1511    
1512    =head3 validate_fasta_file
1513    
1514        $sequence_type = validate_fasta_file($in_file, $out_file)
1515    
1516    Ensure the given file is in valid fasta format. If $out_file
1517    is given, write the data to $out_file as a normalized fasta file
1518    (with cleaned up line endings, upper case data).
1519    
1520    If successful, returns the string "dna" or "protein".
1521    
1522    Will invoke die() on failure; call inside eval{} to ensure full error catching.
1523    
1524    =cut
1525    
1526    sub validate_fasta_file
1527    {
1528        my($file, $norm) = @_;
1529    
1530        my($input_fh, $clean_fh);
1531        open($input_fh, "<", $file) or die "cannot open $file: $!";
1532    
1533        if ($norm)
1534        {
1535            open($clean_fh, ">", $norm) or die "cannot write normalized file $norm: $!";
1536        }
1537    
1538        my $state = 'expect_header';
1539        my $cur_id;
1540        my $dna_chars;
1541        my $prot_chars;
1542    
1543        while (<$input_fh>)
1544        {
1545            chomp;
1546    
1547            if ($state eq 'expect_header')
1548            {
1549                if (/^>(\S+)/)
1550                {
1551                    $cur_id = $1;
1552                    $state = 'expect_data';
1553                    print $clean_fh ">$cur_id\n" if $clean_fh;
1554                    next;
1555                }
1556                else
1557                {
1558                    die "Invalid fasta: Expected header at line $.\n";
1559                }
1560            }
1561            elsif ($state eq 'expect_data')
1562            {
1563                if (/^>(\S+)/)
1564                {
1565                    $cur_id = $1;
1566                    $state = 'expect_data';
1567                    print $clean_fh ">$cur_id\n" if $clean_fh;
1568                    next;
1569                }
1570                elsif (/^([acgtumrwsykbdhvn]*)\s*$/i)
1571                {
1572                    print $clean_fh uc($1) . "\n" if $clean_fh;
1573                    $dna_chars += length($1);
1574                    next;
1575                }
1576                elsif (/^([*abcdefghijklmnopqrstuvwxyz]*)\s*$/i)
1577                {
1578                    print $clean_fh uc($1) . "\n" if $clean_fh;
1579                    $prot_chars += length($1);
1580                    next;
1581                }
1582                else
1583                {
1584                    my $str = $_;
1585                    if (length($_) > 100)
1586                    {
1587                        $str = substr($_, 0, 50) . " [...] " . substr($_, -50);
1588                    }
1589                    die "Invalid fasta: Bad data at line $.\n$str\n";
1590                }
1591            }
1592            else
1593            {
1594                die "Internal error: invalid state $state\n";
1595            }
1596        }
1597        close($input_fh);
1598        close($clean_fh) if $clean_fh;
1599    
1600        my $what = ($prot_chars > 0) ? "protein" : "dna";
1601    
1602        return $what;
1603    }
1604    
1605  sub strip_func {  sub strip_func {
1606          my($func) = @_;          my($func) = @_;
1607    

Legend:
Removed from v.1.38  
changed lines
  Added in v.1.39

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3