LoginSignup
1
0

More than 3 years have passed since last update.

コロナウィルスをネタにBioPerlなプログラミングについて考えてみる

Last updated at Posted at 2020-03-31

はじめに

昨今はネットを通じて多くの学術データが自由に取得可能となり、PCの高性能化やディスクの大容量化も相まって、家庭にあってボランティアベースでできる研究の幅もずいぶん広がってきた。いわゆる「象牙の塔」の惨憺たるありさまとは対照的に、おそらくその辺に従来にはなかったチャンスが潜んでいる。

おりしも、コロナウィルスによる厄災がここ3か月ほどで爆発的に拡大して世界中が大騒動になっており、2020年3月末現在ではそれを収束させるための抜本的な対策がまだ打ち出せていない状況にある。なかなか目が出ず燻っているいるポスドク諸君、バイトすらままならず引きこもっている諸君、今こそネットに、コンピュータに向き合おうじゃないか! 相同性検索だ! 統計処理だ! データマイニングだ! 感染抑止や、もしかしたら創薬なんぞにつながる何らかの知見が我々のPCから見つかるかもしれない!!

…というような気合とはあまり関係なくQiitaをサーフィンしていたところ、タイムリーな話題を取り上げている記事があった。

参照記事:
コロナウィルス の配列を比較するために

簡単なperlスクリプトでとりあえずウィルスの塩基配列データを取ってこようという趣旨の記事である。このような情報共有はますます積極的にすすめるべきであり、まずは記事をアップしてくださっている@ yohoさんに敬意を表したい。

しかしながら、そのスクリプトの中身については少々コメントしたいことがある。本記事ではそれについて述べ、私なりにプログラミングというものについて思うところなどを書いてみたいと思う。

私のコード

細かい話は後回しにして、まずは当該記事に触発されて私が書いてみたperlスクリプトを提示する。当該記事に紹介されていたスクリプト(以下、「元スクリプト」と表記)と(ほぼ)同じ機能を実現するものだ。

#!/usr/bin/env perl

use strict;
use warnings;
use Bio::Seq;
use Bio::SeqIO;

my $infile = 'test.gb';
my $in     = Bio::SeqIO->new(-file => $infile, -format => 'Genbank');

while( my $record = $in->next_seq()){ # Do something for each Genbank record

  # dnaSequence
  my $seq = $record->seq;
  $seq=~/n/i and next;

  # DEFINITION
  my $definition = $record->desc();
  $definition=~/partial/ and next;

  # Output various things
  print "Locus: ",   $record->id,                       "\n";
  print "Species: ", $record->species->scientific_name, "\n";
  foreach my $reference ( $record->annotation()->get_Annotations('reference')){
    ($reference->title())    and print "Title: ",    $reference->title(),    "\n";
    ($reference->authors())  and print "Authors: ",  $reference->authors(),  "\n";
    ($reference->location()) and print "Location: ", $reference->location(), "\n";
  }
  print "Sequcence: $seq\n\n";
}

このスクリプトは、GenBank形式と言われる塩基配列および関連情報を記したファイルから幾種類かのデータを抜き出す。そこは元スクリプトと同じだが、出力は「複行レコード1」の形で出力することにしたほか、出力する内容は元スクリプトに準拠してはいるが微妙に違う。

それ以外の重要な相違点は次のように纏められるだろう。

  • Bioperlのモジュールを利用している。
  • use strict; use warnings; およびmy指定を怠らない。
  • ループの中身は簡素にする。

以下、個々のコードについて、元スクリプトとの比較を交えながら説明し、私の考えを述べることにしよう。

bioperlモジュール

ありもののコードはありがたく使う

何といっても本スクリプトの肝は、bioperlとくにBio::SeqIOとBio::Seqを利用している点にある。これによりスクリプトの内容は元スクリプトと比べて大幅にシンプルなものとなっている。

元スクリプトでは、複数行にまたがるテキストをまとめたりインデントを取ったりする目的で、フラグ変数を使ったり正規表現で加工したりしている。苦労のあとがうかがえる。だが、正直いって見辛い。

腕試しや勉強としてならともかく、実戦ではこんなのは絶対レディメイドの関数に任せてしまったほうが効率的である。

こんなところでエンバグして的外れなデータを次ステップに渡すような事態を招いては目も当てられない。そうしないためにデバッグしようとするとどうしても慎重さが求められるし、時間の消費は避けられず、なにより心の健康に悪い。

モジュールを活用すれば機能拡張も楽である。たとえば、本スクリプトでは生物種を出力する処理を追加した。Bio::Taxonオブジェクトのメソッドを利用できたのでたった1行で片付いている。素のperlコードで対処するとなるとコーディングの手間は大きく増えるはずだ。

ただし「素の」コードが有効なこともある

ただ、私の経験上、文字列データの加工(要するに塩基配列の切り出しやら相補鎖配列の生成など)だけは、bioperlの関数よりは素の(=標準関数のみを用いた)perlコードで対処したほうが早い。標準関数を使って一行で書けるような処理はそのように書こう。無理やりオブジェクト指向なカッコいいが回りくどいコードを書く必要はない。

オブジェクトやメソッドに関するあれこれ

Bio::SeqIOはGenBank形式だけでなくFastaをはじめ様々な形式に対応しており、どの形式に対しても同じような手順で必要なデータにアクセスできるように設計されている。そうした汎用性と引き換えに、やや自分の直感に反する在り様のメソッドを利用しなければならないこともでてくる。たとえばGenBank形式におけるDEFINITION行を持ってくるのがdefinitionではなくdescだったりして、なんか落ち着かなく感じるのだが、こういうケースではもうperldoc Bio::SeqIOで説明をチェックしまくり「そういうもんだ」と何とか納得する以外ない。

このように、モジュールは便利だが、特にオブジェクト指向バリバリなモジュールの場合、作者の世界観にどっぷりつかり受け入れなければコードを使いこなせないきらいがある。その点がちょっとストレスフルだったりはする。それを嫌ってモジュールを使わない方向に行くのも一つの方法ではあるし、前述のようにそれが有効な場合もあるのは確かだ。

しかし、ある程度以上複雑な処理なら、やはり普通はモジュール作者の軍門に下ってその世界観を受け入れ、ありがたく使わせてもらったほうがいい。

use strict・use warnings・my指定

使い捨てのワンライナーならともかく、仕事でがっつり使おうというしっかりしたスクリプトを作りたいなら、冒頭に

use strict;
use warnings;

を書くのは、マストである。その上で、すべての変数はmy $x;のようにmyをつけた使用宣言を行う2。一部界隈から不満が聞こえそうだ。「カジュアルに自由に変数が使えるからperlなんだよ。面倒じゃん。いらねーよ。」
そのような主張には全力で否といわせてもらう。繰り返すがこの宣言はマストである。なぜなら、それだけで変数名のスペルミスに起因する間抜けなエンバグを防げるのだ。

この宣言をしないでいきなり変数を使おうとした場合、perlはその変数の中身を適当に0とか''とかだと仮定して処理を進めようとする。ある場所で$seq="ACAC"と書くべきところ、間違えて$req="ACAC"と書いてしまったとしよう。以降で$seqを使った全ての処理が間違った値を吐き出すことになる。だがperlはその場ではエラーを出さず、無意味な処理を黙って続けようとするはずだ。ずっと先で破綻をきたして止まるか、あるいは的外れな答えをもっともらしく出力して止まる。この状況がどれほど疎ましいかは説明しなくてもいいだろう。myを使うだけでそれを排除できるのである。使わない手があるだろうか?

ループの中身はできるだけ簡素にする

forとかwhileの下に続くコードはできるだけ短くしよう。

とくに多重ループを使っている場合、一番内側のループで無駄なことをしないよう殊更しっかりリファクタリングするべきだ。ループの外でできること(定数の定義など)は全部ループの外に出す。それだけで全体の所要時間があっというまに1/10とかになったりするぞ。

本スクリプトの中で... and nextというコードを使っているのを見てほしい。元スクリプトでは「塩基配列中に文字nが含まれている」場合、そのデータは棄却して処理を行わないことになっている。そのチェックのために「フラグを返す関数」を利用しているのだが、こういうケースではメインループの直下に条件判定を書いてしまうのが手っ取り早い。

排除条件式 and next;

という形で、「条件が成立した場合はその場で当該データの処理は中止して次へ」というコードを入れる。nextが実行されると、その場で実行は中断されループの先頭に制御が戻って、ループの次の回の処理に移行する。場合によってはnextの代わりにlastが使えるだろう。その場合ループの周回自体が中止され、ループの後の処理に移る。

このようなコードをループ内の可能な限り上の方に書くのが望ましい。必要ないとわかったループからは可及的速やかに離脱するのだ。

もちろんifブロックを使ってもいい。しかし、nextを使ったほうがシンプルで、ネストを深くせずに済み、大いにコードの保守性が上がる。

おわりに

誤解してほしくないのだが、私に元記事の価値を否定する意思は毛頭ない。

ただ、ちょっともったいなくも思ったので、問題提起のための便乗記事を書かせてもらった。もしよければコーディングの心得について私なりの見解を述べた過去記事『perl: system関数への丸投げはヤメよう』も見ていただけると嬉しい。


  1. 空行で区切られたブロックを一つのレコードと見なし改行をフィールドの区切りと見なす形式。 

  2. 厳密にはourというのもあるのだが、まあ使う機会はほとんどないだろう。原則myと覚えておいて差し支えない。 

1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
1
0