RやPythonであれば、さまざまな洗練された既存のライブラリの使い方について記載することになると思うけれども、Rubyでは既存のライブラリが多くないので、もっぱら自分でC言語のライブラリのバインディングを作成し、それを使う感じになる。
2bitファイルを使う場合
UCSCのサイトから.2bitファイル形式のゲノム配列を入手することができる。たとえば、
https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/
から hg38.2bit をダウンロードすることができる。
bio-twobit というライブラリを作成したのでそれを使う。これはC言語のライブラリ lib2bit のバインディングであり高速に動作する。
require 'bio/twobit'
hg38 = Bio::TwoBit.open("hg38.2bit")
hg38.sequence("chr1", 50000, 50050)
# "AAACAGGTTAATCGCCACGACATAGTAGTATTTAGAGTTACTAGTAAGCC"
hg38.close
ここで、注意するべきなのは、hg38.sequence("chr1", 50000, 50050) の最初の引数は0-baseで、2番目の引数は1-baseであることである。あるいは、左閉右開の半開区間と考えることもできる。この仕様は、元プロジェクトのpy2bitに準じている。
また、このツールは、多くのリファレンスゲノムに対応しており、
hg19 = Bio::TwoBit::Hg19.new
hg38 = Bio::TwoBit::Hg38.new
hs1 = Bio::TwoBit::Hs1.new
などとすると、自動的に2bit形式のゲノムをダウンロードして、呼び出せるようになります。ダウンロードしたファイルはキャッシングされます。
Fastaファイルを使う場合
ruby-htslib というライブラリを作成したのでそれを使う。
これはC言語でHTSファイルを扱うための標準的なライブラリであるhtslibのRubyバインディングである。
require "htslib"
f = HTS::Faidx.new("GRCh38.p13.genome.fa")
f.seq("chr1", 50000, 50049)
# => "AAACAGGTTAATCGCCACGACATAGTAGTATTTAGAGTTACTAGTAAGCC"
f.close
ここでは、f.seq("chr1", 50000, 50049) において、いずれの引数も0-baseであることに注目してほしい。
これはhtslibのネイティブ関数 faidx_fetch_seq の仕様に従っている。
注意してもらいたいのは、chr1:50001-50050 のように範囲を文字列で表すときは、1-based の閉区間になることだ。
require "htslib"
f = HTS::Faidx.new("GRCh38.p13.genome.fa")
f.seq("chr1:50001-50050")
# => "AAACAGGTTAATCGCCACGACATAGTAGTATTTAGAGTTACTAGTAAGCC"
f.close
Samtoolsのfaidxコマンドでも
samtools faidx GRCh38.p13.genome.fa "chr1:50001-50050"
# >chr1:50001-50050
# AAACAGGTTAATCGCCACGACATAGTAGTATTTAGAGTTACTAGTAAGCC
とすると、同様のじ結果が表示される。もちろん、これを open3 や system などで呼び出してもいい。というか普通はこの方法が第一選択だと思われる。
とにかく、0-based と 1-based は難しい。
すぐに間違える。迷ったらその都度確認する方がよさそうである。
この記事は以上です。