Getting the wt/mutant AA based on genomic SNP position

get_aa.rb

require 'rubygems'
require '/nfs/team29/ruby-ensembl-api/lib/ensembl'

codon_table = Bio::CodonTable[1]

include Ensembl::Core
DBConnection.connect('homo_sapiens', 52)

STDIN.each do |line|
  gene_name, transcript_acc, chr, pos, ref, alleles = line.chomp.split(/\t/)
  transcript = Transcript.find_by_stable_id(transcript_acc)
  ref = ref.downcase
  alleles = alleles.downcase
  alt_allele = alleles.sub(/\//,'').sub(/#{ref}/,'').split(//)[0]
  if transcript.strand == -1
    ref.tr!('acgt','tgca')
    alt_allele.tr!('acgt','tgca')
  end
  begin
    cds_pos = transcript.genomic2cds(pos.to_i)
  rescue
    STDERR.puts transcript_acc + "\t" + $!
    next
  end
  codon_nr = (cds_pos.to_f/3).ceil
  frame = cds_pos%3 - 1
  frame = 3 if frame == -3
  start_of_codon = 3*codon_nr - 2
  codon = transcript.cds_seq.split(//)[start_of_codon-1,3].join('')
  mutated_codon = codon.clone
  mutated_codon[frame] = alt_allele
  STDERR.puts "WARNING: wrong AA" unless transcript.protein_seq.split(//)[codon_nr-1] == codon_table[codon]
  puts [chr + '_' + pos, 'undef', transcript_acc, codon_nr, codon_table[codon], codon_table[mutated_codon]].join("\t")
end