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ゼミ(サンプルコード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しているかどうか比較する、簡単なコードを書くことになるのであろう。