utils/human_mrna_cds_insert.pl

   1 #!/usr/bin/perl
   2
   3 use warnings;
   4 use strict;
   5
   6 # psql snp -c "COPY mrna_cds_table FROM STDIN WITH DELIMITER '   ' NULL AS 'NULL'";
   7
   8 my $prot_table;
   9
  10 while (<DATA>) {
  11     chomp;
  12     my ($codon,$aa,$aa_med,$aa_long) = split /\t/;
  13     $prot_table->{$codon} = $aa;
  14 }
  15
  16 my $current_stanza;
  17 while (<>) {
  18     $current_stanza .= $_;
  19     if ($_ =~ m{^//$}) {
  20         handle_mrna_stanza($current_stanza);
  21         $current_stanza = '';
  22     }
  23 }
  24
  25 sub handle_mrna_stanza{
  26     my ($current_stanza) = @_;
  27
  28     my ($accession,$version,$gi) =
  29         $current_stanza =~ m/^VERSION\s+([^\.]+)\.(\d+)\s+GI:(\d+)$/m;
  30     my ($gene) =
  31         $current_stanza =~ m/^\s+\/gene=\"?([^"]+)\"?$/m;
  32     my ($cds) =
  33         $current_stanza =~ m/^\s+CDS\s+(\d+(?:\.\.\d+))$/m;
  34     return if not defined $cds;
  35     my ($start,$stop) = split /\.\./,$cds;
  36     return if not defined $start;
  37     if (not defined $stop) {
  38         $stop = $start;
  39     }
  40     my ($sequence) =
  41         $current_stanza =~ m/^(ORIGIN.+)$/sm;
  42     $sequence =~ s{\s*//\s*$}{};
  43     $sequence =~ s{^ORIGIN}{};
  44     $sequence =~ s{[\d\s\n]+}{}g;
  45     $sequence = uc($sequence);
  46     my $sequence2 = substr $sequence,($start-1),($stop-$start+1);
  47     # There are other posibilities, like CTG, ACG, etc; so we can't test this.
  48     # if ($sequence2 !~ /^A[UT]G/i) {
  49     #   print STDERR $sequence2."\n";
  50     #   print STDERR $gene."\n";
  51     #   die "Sequence doesn't start with a start codon";
  52     # }
  53     if ($sequence2 !~ /(?:[TU]AA|[TU]AG|[TU]GA)$/) {
  54         print STDERR $sequence2."\n";
  55         print STDERR $gene."\n";
  56         die "Sequence doesn't end with a stop codon";
  57     }
  58     # figure out cds gi
  59     my $cds_gi='NULL';
  60     my ($cds_stanza) = $current_stanza =~ /^\s{5}CDS\s+(.+?)^(?:\s{5}\w+|ORIGIN)/ms;
  61     if (length $cds_stanza) {
  62         ($cds_gi) = $cds_stanza =~ /db_xref="GI:(\d+)"/;
  63         if (not defined $cds_gi) {
  64             $cds_gi = 'NULL';
  65         }
  66     }
  67     print join("\t",$gi,$cds_gi,$accession,$version,$gene,$start,$stop,length($sequence2), uc($sequence2),convert_to_prot_seq($sequence2))."\n";
  68 }
  69
  70 sub convert_to_prot_seq{
  71     my ($mrna) = @_;
  72 $mrna =~ s/U/T/g;
  73     my $prot = '';
  74     for my $pos (0..(length($mrna) / 3 - 1)) {
  75         my $codon = substr($mrna,$pos*3,3);
  76         if (not exists $prot_table->{$codon}) {
  77            $prot .= 'N';
  78         }
  79         else {
  80             $prot .= $prot_table->{$codon};
  81         }
  82     }
  83     return $prot;
  84 }
  85
  86 __DATA__
  87 ATT     I       Ile     Isoleucine
  88 ATC     I       Ile     Isoleucine
  89 ATA     I       Ile     Isoleucine
  90 CTT     L       Leu     Leucine
  91 CTC     L       Leu     Leucine
  92 CTA     L       Leu     Leucine
  93 CTG     L       Leu     Leucine
  94 TTA     L       Leu     Leucine
  95 TTG     L       Leu     Leucine
  96 GTT     V       Val     Valine
  97 GTC     V       Val     Valine
  98 GTA     V       Val     Valine
  99 GTG     V       Val     Valine
 100 TTT     F       Phe     Phenylalanine
 101 TTC     F       Phe     Phenylalanine
 102 ATG     M       Met     Methionine
 103 TGT     C       Cys     Cysteine
 104 TGC     C       Cys     Cysteine
 105 GCT     A       Ala     Alanine
 106 GCC     A       Ala     Alanine
 107 GCA     A       Ala     Alanine
 108 GCG     A       Ala     Alanine
 109 GGT     G       Gly     Glycine
 110 GGC     G       Gly     Glycine
 111 GGA     G       Gly     Glycine
 112 GGG     G       Gly     Glycine
 113 CCT     P       Pro     Proline
 114 CCC     P       Pro     Proline
 115 CCA     P       Pro     Proline
 116 CCG     P       Pro     Proline
 117 ACT     T       Thr     Threonine
 118 ACC     T       Thr     Threonine
 119 ACA     T       Thr     Threonine
 120 ACG     T       Thr     Threonine
 121 TCT     S       Ser     Serine
 122 TCC     S       Ser     Serine
 123 TCA     S       Ser     Serine
 124 TCG     S       Ser     Serine
 125 AGT     S       Ser     Serine
 126 AGC     S       Ser     Serine
 127 TAT     Y       Tyr     Tyrosine
 128 TAC     Y       Tyr     Tyrosine
 129 TGG     W       Trp     Tryptophan
 130 CAA     Q       Gln     Glutamine
 131 CAG     Q       Gln     Glutamine
 132 AAT     N       Asn     Asparagine
 133 AAC     N       Asn     Asparagine
 134 CAT     H       His     Histidine
 135 CAC     H       His     Histidine
 136 GAA     E       Glu     Glutamic acid
 137 GAG     E       Glu     Glutamic acid
 138 GAT     D       Asp     Aspartic acid
 139 GAC     D       Asp     Aspartic acid
 140 AAA     K       Lys     Lysine
 141 AAG     K       Lys     Lysine
 142 CGT     R       Arg     Arginine
 143 CGC     R       Arg     Arginine
 144 CGA     R       Arg     Arginine
 145 CGG     R       Arg     Arginine
 146 AGA     R       Arg     Arginine
 147 AGG     R       Arg     Arginine
 148 TAA     *       *       Stop codon
 149 TAG     *       *       Stop codon
 150 TGA     *       *       Stop codon