Perlの参照を使いBioMart出力をまとめる

発現変動のみられた遺伝子群のリストでアノテーションをBioMartから抽出してみた。
するとこういう状態になってしまっている。

id	scheme1	scheme2
gene000001	DescAA	DescBA
gene000001	DescAA	DescBB
gene000001	DescAA	DescBC
gene000001	DescAB	DescBA
gene000001	DescAB	DescBB
gene000001	DescAB	DescBC
gene000002	DescAC	DescBD
gene000002	DescAC	DescBE
gene000002	DescAC	DescBF
gene000002	DescAD	DescBD
gene000002	DescAD	DescBE
gene000002	DescAD	DescBF

これを

id	scheme1	scheme2
gene000001	DescAA, DescAB	DescBA, DescBB, DescBC
gene000002	DescAC, DescAD	DescBD, DescBE, DescBF

のようになおしたいのだった。

そんなことを後ろに座っている院生と話すと、彼は、自分はRで書いてみるから競争しよう、ということで彼はR、私はPerlを使って試みるプチハッカソンが午後突如始まった。

プログラミングは例によって不調法なのだが、Perlの配列やハッシュの「参照」というのを使えばできるらしいという考えはあった。

私はPerlを結城浩先生の名著Perl言語プログラミングレッスンで学んだ。

新版Perl言語プログラミングレッスン入門編

この本がなければ自分はPerlどころかプログラムを書くという構えにすら到達できなかったのは間違いない。考え方そのものを学んだ大切な本だった。

ただ、私はこのとき、参照というものを学んでいない。

マインドマップもかいてかなり読み込んだと思っていた。そのときのマインドマップがこれである。クリックすると元画像が開くと思う。

一念発起して参照に関する記述と首っ引きで作った。
参考 リファレンスの使い方をマスターしよう - Perlゼミ(サンプルコードPerl入門)

もちろん、結城先生ので勉強しなければこうしてつまみ読みしてつくることすらおぼつかなかったので、感謝してもしきれません。

で、まず出力側から試してみた。

#!/usr/bin/perl
use strict;
use warnings;

# 出力側コード

# 各遺伝子をハッシュのキーに、そのハッシュの内容は配列になっていて、重複を含んだハッシュキー群をリストしている
my $ids = {
	'gene000001' =>
		[ { 'DescAA' => 1, 'DescAB' => 2, 'DescAB' => 3 },
		{ 'DescBA' => 1, 'DescBA' => 2, 'DescBB' => 3, 'DescBC' => 4 }
		],
	'gene000002' => [
		{ 'DescAC' => 1, 'DescAC' => 2, 'DescAD' => 3 },
		 { 'DescBD' => 1, 'DescBE' => 2, 'DescBF' => 3, 'DescBF' => 4 }
		 ],
	};


for my $id (sort keys %$ids ) { # geneごとに取り出す
	print "$id\t"; # gene名前を表示してタブをおく
	
	my $schemecount = @{$ids->{$id}}; # いくつ枠 (GO, KOG, SMART, etc.) があるか
	for (my $i = 0; $i < $schemecount; $i++) { # ひとつずつ取り出す
		my @scheme =  keys %{$ids->{$id}->[$i]}; # 各枠に入ったハッシュのキーを配列に入力
		print join(' | ', @scheme); # バーで区切って表示
		print "\t"; # タブで枠間を区切る
	};
	print "\n"; # 改行する
}

で、これでいちおううまくいくのが確認できたので、入力側も作ってみた。

#!/usr/bin/perl
use strict;
use warnings;

my $file = "2015MAR24.txt";

open (IN, $file) or die "$!";


my $genes = {};
while (my $line = <IN>) {
	chomp $line;
	
	my ($gene, @annotation) = split('\t', $line);

	$genes->{$gene} ||= [];
	
	my $annotationnumber = @annotation;
	for (my $i = 0; $i < $annotationnumber; $i++) {
		$genes->{$gene}->[$i] ||= {};
		
		$genes->{$gene}->[$i]{$annotation[$i]} = 1;
	}
};

for my $id (sort keys %$genes ) {
	print "$id\t";
	
	my $schemecount = @{$genes->{$id}};
	for (my $i = 0; $i < $schemecount; $i++) {
		delete $genes->{$id}->[$i]{""}; # 空白キーの削除
		my @scheme =  keys %{$genes->{$id}->[$i]};
		print join(' | ', @scheme);
		print "\t";
	};
	print "\n";
}

これにて無事可能になった。

今回は記述をまとめるコードを作ったことになる。次は、functional annotation IDの数えてenrichmentしているかどうか比較する、簡単なコードを書くことになるのであろう。