遺伝子のプロモータ配列を簡単に得るスクリプト

遺伝子のプロモータを得るだけのために、各種ゲノム配列を準備するのは面倒。NCBI のE-Utils利用すると、ゲノム配列の準備とかしなくても、プロモータ配列が得られます*1。

準備するもの

• ruby
• rubygems
  • xmlsimple
• ネットワーク環境
• プロモータを得たい遺伝子のNCBI Gene ID(群)

使い方
引数にプロモータを得たい遺伝子のNCBI Gene IDをカンマ区切りで列挙します。
以下の様にすると、3つのIgG (FCGR1A(2212), FCGR2B(2213),FCGR3A(2214))のプロモータを取得して表示します(転写開始点から上流250bp, 転写開始点下流50bpの計300bp)。

% ruby get_promoter_seqs.rb 2212,2213,2214

取る領域を変更したい場合はスクリプト内のPROMLEN(取る配列の長さ)とPROMINGENE(転写開始点から下流を何塩基取るか)を変更してください。

結果のサンプル
MultiFasta形式で、注釈にはNCBI Gene ID, Symbol名, 取得した配列のゲノム上の場所が書かれます。

% ruby get_promoter_seqs.rb 2212,2213,2214
>2212|FCGR2A| from >ref|NC_000001.9|NC_000001:159741593-159741892 Homo sapiens chromosome 1, reference assembly, complete sequence
CCGCGCCCGGCCGACCAGGAATTATCTTAAACAGAAAAAGAATACTCTAAGGAGGGGTATACCGGCCAGC
AGAACAGTAACCCCTCCCCGGGGTTTTCACAAGGACACAAGCTGGTGTGTAATAGACTTATCTATAAAGC
AGACTCAACTCCTAACTGGAATTGTTCCTGATTGGTTAGCAGAAACAGAATTGAGGACTGACGACAGCTG
CACAAGAAGATGAAACCATTTCCTTCCTCTTTTCTAAGCTTGTCTCTTAAAACCCACTGGACGTTGGCAC
AGTGCTGGGATGACTATGGA
>2213|FCGR2B| from >ref|NC_000001.9|NC_000001:159899313-159899612 Homo sapiens chromosome 1, reference assembly, complete sequence
CTGAGCTGAGTTTTGGTGACCATTTCCCTCTTTCTCCCAGAGGCCCAGGCCAGCTGTGGCCTCAGAGGAA
GAAGAAGGGAGTTGTTTCCCTAGTTTCTAAAATTTCTGTGAATTTGAACATGGGCTACACCAGATTTATT
CTGGGAAGCTCTGAATCTTCTAGGAGGGAAAGACTGAGAGGAAAGAGGGTGGAAAGGGAGGAGCCTGTGA
TAAAACAGAACATTTCTTTTTCACTTCCCCTTTCAGACTCCAGAATTTGTTTGCCCTCTAGGGTAGAATC
CGCCAAGCTTTGAGAGAAGG
>2214|FCGR3A| from >ref|NC_000001.9|NC_000001:159786392-159786691 Homo sapiens chromosome 1, reference assembly, complete sequence
TTTACTCCCTCAAAGGTCTGTGGCTGAGCATCTGAGGACACACACAGAATCTGCCAGAGTGTGCCCTCAG
CCATCCCAGGATGCTTGCCCCATCTCCTGGATTTGGATCCACCCAGCACCAAGAATGGGACAGTGAGACC
CTGGGGATGAGATTCAAGGTGGGAGGAGCATTCTCTGAGGGCTTCCTTCATTTCACCTCAGAACTTGTCA
CTCTCCTGCCTCGCCCAGACCCATCTATCTCCAGCTGAGGCCCTGCCTGCACACAGAAAGTGGTCCTTTC
AAATCTTCAGACCACTAGCA

スクリプト
get_promoter_seqs.rb

#!/usr/bin/env ruby

require 'rubygems'
require 'xmlsimple'
require 'open-uri'
require 'pp'

PROMLEN = 300 # length of promoter
PROMINGENE = 50

# argument is a list of Entrez gene ids separated by commas
gids = ARGV.shift

# read gene information summary
uri = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?"
uri += "db=gene&id=#{gids}&retmode=xml"
#uri = "gid_780_5982.xml"
xmldata = open(uri).read
genesinfo = XmlSimple.xml_in(xmldata)

# parse gene information xml
genesinfo["DocSum"].each do |ginfo|
  chrstart = nil
  chrstop = nil
  chraccver = nil
  geneid = ginfo["Id"][0].to_i
  genename = nil
  begin
    ginfo["Item"].each do |item|
      if item["Name"] == "Name"
        genename = item["content"]
      elsif item["Name"] == "GenomicInfo"
        chrinfo = item["Item"][0]["Item"]
        chrinfo.each do |info|
          if info["Name"] == "ChrStart"
            chrstart = info["content"].to_i
          elsif info["Name"] == "ChrStop"
            chrstop = info["content"].to_i
          elsif info["Name"] == "ChrAccVer"
            chraccver = info["content"]
          end
        end
      end
    end
  rescue
    STDERR.puts "Read Error at #{geneid}"
    next
  end
  promstart = nil
  promstop = nil
  strand = nil
  if chrstart < chrstop 
    promstart = chrstart + PROMINGENE - PROMLEN
    promstop = chrstart + PROMINGENE - 1
    strand = 1
  else
    promstart = chrstart - PROMINGENE + PROMLEN
    promstop = chrstart - PROMINGENE + 1
    strand = 2
  end
  
  query = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?"
  query += "db=genome&id=#{chraccver}&seq_start=#{promstart}"
  query += "&seq_stop=#{promstop}&strand=1" # "#{strand}"
  query += "&rettype=fasta"

  # get its promoter sequence
  fp = open(query)
  promseq = fp.read().split("\n")
  fp.close
  # display the promoter sequence
  promseq[0] = ">" + geneid.to_s + "|" + genename.to_s + "| from " + promseq[0] 
  puts promseq.join("\n")
end

*1:TogoWS使いたいんだけど、部分配列の抽出に対応していないので断念対応しているそうです。コメント参照。