RSS「オープンソース」を使ってみよう (第27回 フリーソフトウェアを用いたゲノム科学におけるビッグデータ処理)
01/24
Use it ! OSS No Comments
Tweet
—————————————————————————-
東京農工大学ゲノム科学人材育成プログラム特任教授
石井 一夫
—————————————————————————-
1. ゲノム科学で用いられるフリーソフトウェア
次世代シーケンサーというDNA塩基配列情報を大量に産生する機器が実用化
されて数年が経過し、ゲノム科学を扱う医学、生物学の世界では日々のDNA
塩基配列データの産生量や、その取り扱うデータ量が飛躍的に増えています。
これらのデータ処理にはUNIX/Linuxを中心とするフリーソフトウェアは欠か
せません。
(1) 汎用のフリーソフトウェア
次世代シーケンサーのデータは、例えばイルミナ社製の解析機器から産生
されるデータの場合、1ファイルあたりに数千万断片から数億断片のDNA
塩基配列データとそのクウォリティデータを含むファイルが産生されます。
それを、(1) catやgrep、sed、awkなどのシェルのコマンドや、(2) Perl、
Python、Rubyなどのスクリプト言語、(3) R、Octaveなどの統計解析言語
を組み合わせて処理します。
必要に応じて、(4) MySQL、PostgreSQLなどのデータベースも使用します。
むしろ、このようなスケールのデータ解析では、データベースは必須です。
(2) 生物学的なデータ解析専用のソフトウェア
もちろん、生物学的な情報解析専用のソフトウェアも多数開発されています。
例えば、(1) 次世代シーケンサーから産生されたDNA塩基配列データを互いに
ジグソーパズルのように結合させ、長いDNA塩基配列を得るアセンブリと呼ば
れるデータ解析行程では、Velvet、Oases、Trinityなどが使用されます。
また、(2) 次世代シーケンサーから産生されたDNA塩基配列データを既知の
DNA塩基配列へ整列させるマッピングと呼ばれるデータ解析行程では、BWA、
Bowtieなどが使用されます。
(3)このようにして得られた長めのDNA塩基配列データと既知のDNA塩基配列
との相同性の検索には、BLASTなどが用いられます。Perl、Python、Ruby、
Javaなどには、(4) 生物学的なデータの解析に特化した関数などを集めた
ライブラリが整備されており、それぞれBioPerl、BioPython、BioRuby、
BioJavaと呼ばれています。
また、統計解析ソフトのRには、(5) Bioconductorと呼ばれる生物学的解析用
のパッケージ群が存在します。
(3) ビッグデータ処理に特化したソフトウェア
近年、(1) クラウドコンピューティングや仮想化環境を組み合わせたビッグ
データ処理の技術がゲノム科学にも応用されており、Hadoop、OpenStack、
Eucalyptusなどが使われています。ゲノム科学においてはビッグデータ処理
のためのファイルシステムも注目されており、(2) FreeBSDにおけるZFSなど
ペタバイト級、エクサバイト級のスケールのビッグデータ処理も視野に入れ
た次世代ファイルシステムの導入も進んでいます。私たちの教育プログラム
でゲノムデータ解析に使用している主なフリーソフトを表１に示します。
表：東京農工大学ゲノム科学人材育成プログラムでデータ解析に用いるソフトウェア群
Table：Software for Genomic Data Analysis used in Tokyo University of Agriculture and Technology
|
2. ゲノム科学におけるビッグデータ解析
現在は、ゲノム科学で扱うDNA塩基配列データのファイルのサイズは数百メガ
バイトから数十ギガバイトのサイズのものが多く、これはファイルの行数で表現
すると数万行から数億行くらいです。
これらをPerl、Python、Rubyなどのスクリプト言語やR、Octaveなどの統計解析
ソフト、BLASTなどのDNA塩基配列解析の専用ソフトウェアでデータ処理をする
場合に、メモリ不足となりソフトウェアの実行動作が極端に遅くなったり、ハン
グアップを起こしたりということが頻繁に発生します。
メモリを増やすなどのハード面の充実によりこれらの解析を可能にすることも
できますが、ファイルを分割して単位プロセス当たりのデータ量を小さくし、
並列処理やバッチ処理による直列処理を組み合わせることにより、効率的に解
析を行えます。近年は、コア数が多いCPUを搭載したコンピュータが増えてい
るので、それを活用した解析戦略が有効だと思われます。
3. ゲノム科学におけるビッグデータ解析の例
ゲノムデータでよく行われる作業は、ファイル形式の変換です。
ここでは、動作はそれほど速くないので、あまり実用的ではありませんが、
ひとつの練習例としてRとPerlをシェルスクリプトと組み合わせてGFFという
ファイル形式から、遺伝子の塩基配列データを収納したマルチfastaという形式の
ファイルを作成する方法例をご紹介します。
ここでご紹介する方法は最善の例ではありませんし、Perlの代わりにRubyやPythonを
使うことも可能でしょう。
また、ひとつのプログラミング言語ですべてのプロセス完結させることも可能
ですし、同様の操作を行ってくれるツールキット（既存のソフトウェアパッケージ）
も存在すると思います。
しかし、そこまでやらなくても、日常の何気ない作業のなかで、特に特殊なアプリ
ケーションを使用しなくても短時間でさらりと処理できて便利です。
(1) ゲノムDNAの塩基配列情報とアノテーション情報の入手
ここでは、モデル生物であるシロイヌナズナのゲノムDNA塩基配列データと
（遺伝子の位置情報を示す）アノテーションデータから、各遺伝子のマルチ
fasta形式のファイルを作成します。
作成に必要なデータを含むファイルを以下のようにして入手します。
シロイヌナズナのゲノムDNA塩基配列データの入手
$ wget ftp://ftp.arabidopsis.org/home/tair/Sequences/whole_chromosomes/TAIR10*
シロイヌナズナの（遺伝子の位置情報を示す）アノテーションデータの入手
$ wget ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff
以下のように得られたファイルが確認できます。
$ ls TAIR10_GFF3_genes.gff TAIR10_chr2.fas TAIR10_chr4.fas TAIR10_chrC.fas TAIR10_chr1.fas TAIR10_chr3.fas TAIR10_chr5.fas TAIR10_chrM.fas
(2) ゲノムDNA塩基配列情報をデータ抽出しやすいように加工
まず、染色体ごとのfasta形式のファイルをコマンドcatでひとつのファイル
”TAIR10.fas”にまとめます。
$ cat *fas | sed -e '/^$/d' >> TAIR10.fas $ ls TAIR10.fas TAIR10_chr2.fas TAIR10_chr5.fas TAIR10_GFF3_genes.gff TAIR10_chr3.fas TAIR10_chrC.fas TAIR10_chr1.fas TAIR10_chr4.fas TAIR10_chrM.fas
TAIR10.fasを、＜配列名＞タブ＜配列＞のような形に成形するために、
以下のようなPerlスクリプト”getIDseq.pl”を準備しておきます。
$ vi getIDseq.pl
#!/usr/bin/perl
# name: getIDSeq.pl
# usage:
# perl getIDSeq.pl input_filename > output_filename
$input = $ARGV[0];
open ( FILEHANDLE , "< $input");
@array = <FILEHANDLE> ;
chomp(@array);
for (my $i = 0; $i <= $#array; $i++){
if ($array[$i] =~ /^>/) {
if ($i == 0) {
print "$array[$i]\t";
} else {
print "\n$array[$i]\t";
}
} else {
print "$array[$i]";
}
}
print "\n";
close(FILEHANDLE);
以下のコマンドを実行することにより、＜配列名＞タブ＜配列＞の
ような形のファイル“ID_seq.txt”ができます。
$ perl getIDseq.pl TAIR10.fas | awk -F"\t" '{print substr($1,2,4)"\t"$2}' | sed -e 's/chlo/ChrC/g' | sed -e's/mito/ChrM/g' > ID_seq.txt
“ID_seq.txt”は、以下のような形式のファイルになります。
染色体名と塩基配列の間はタブで区切られています。
Chr1 <シロイヌナズナの第1番染色体ゲノムの塩基配列>
Chr2 <シロイヌナズナの第２番染色体ゲノムの塩基配列>
Chr3 <シロイヌナズナの第３番染色体ゲノムの塩基配列>
Chr4 <シロイヌナズナの第４番染色体ゲノムの塩基配列>
Chr5 <シロイヌナズナの第５番染色体ゲノムの塩基配列>
ChrC <シロイヌナズナの葉緑体ゲノムの塩基配列>
ChrM <シロイヌナズナのミトコンドリアゲノムの塩基配列>
(3) アノテーション情報から、転写産物の位置情報を抽出
アノテーションファイル”TAIR10_GFF3_genes.gff”は以下のようになって
います。実のところ、これからコマンドgrepかawkを使って“gene”という
単語を含む行を集めれば簡単に遺伝子の開始位置と終了位置の位置情報を
取り出せます。しかし、今回は練習の意味もあり、exonの情報しかなかった
と仮定して、その情報から各遺伝子の開始位置と終了位置の情報を取り出す
方法をご紹介します。
$ head TAIR10_GFF3_genes.gff Chr1 TAIR10 chromosome 1 30427671 . . . ID=Chr1;Name=Chr1 Chr1 TAIR10 gene 3631 5899 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010 Chr1 TAIR10 mRNA 3631 5899 . + . ID=AT1G01010.1;Parent=AT1G01010;Name=AT1G01010.1;Index=1 Chr1 TAIR10 protein 3760 5630 . + . ID=AT1G01010.1-Protein;Name=AT1G01010.1;Derives_from=AT1G01010.1 Chr1 TAIR10 exon 3631 3913 . + . Parent=AT1G01010.1 Chr1 TAIR10 five_prime_UTR 3631 3759 . + . Parent=AT1G01010.1 Chr1 TAIR10 CDS 3760 3913 . + 0 Parent=AT1G01010.1,AT1G01010.1-Protein; Chr1 TAIR10 exon 3996 4276 . + . Parent=AT1G01010.1 Chr1 TAIR10 CDS 3996 4276 . + 2 Parent=AT1G01010.1,AT1G01010.1-Protein; Chr1 TAIR10 exon 4486 4605 . + . Parent=AT1G01010.1
アノテーションファイル”TAIR10_GFF3_genes.gff”のうちexonの情報だけを
取り出す場合は以下のようなコマンドを実行します。
$ cat TAIR10_GFF3_genes.gff | awk '($3=="exon"){print}' | head
Chr1 TAIR10 exon 3631 3913 . + . Parent=AT1G01010.1
Chr1 TAIR10 exon 3996 4276 . + . Parent=AT1G01010.1
Chr1 TAIR10 exon 4486 4605 . + . Parent=AT1G01010.1
Chr1 TAIR10 exon 4706 5095 . + . Parent=AT1G01010.1
Chr1 TAIR10 exon 5174 5326 . + . Parent=AT1G01010.1
Chr1 TAIR10 exon 5439 5899 . + . Parent=AT1G01010.1
Chr1 TAIR10 exon 8571 8737 . - . Parent=AT1G01020.1
Chr1 TAIR10 exon 8417 8464 . - . Parent=AT1G01020.1
Chr1 TAIR10 exon 8236 8325 . - . Parent=AT1G01020.1
Chr1 TAIR10 exon 7942 7987 . - . Parent=AT1G01020.1
上記のデータのうち今回の解析に必要な情報は、遺伝子名と染色体名、各exonの
開始位置と終了位置です。
その情報を以下のコマンドを用いて抽出し、”exon_list.txt”として保存します。
$ cat TAIR10_GFF3_genes.gff | awk '($3=="exon"){print $9"\t"$1"\t"$4"\t"$5}' | sed -e 's/Parent=//g' | head
AT1G01010.1 Chr1 3631 3913
AT1G01010.1 Chr1 3996 4276
AT1G01010.1 Chr1 4486 4605
AT1G01010.1 Chr1 4706 5095
AT1G01010.1 Chr1 5174 5326
AT1G01010.1 Chr1 5439 5899
AT1G01020.1 Chr1 8571 8737
AT1G01020.1 Chr1 8417 8464
AT1G01020.1 Chr1 8236 8325
AT1G01020.1 Chr1 7942 7987
$ cat TAIR10_GFF3_genes.gff | awk '($3=="exon"){print $9"\t"$1"\t"$4"\t"$5}' | sed -e 's/Parent=//g' > exon_list.txt
”exon_list.txt”から以下のコマンドを用いて、
遺伝子名だけを抽出し、”gene_name.list”という名前で保存します。
$ cat exon_list.txt | cut -f1 | sort | uniq | head AT1G01010.1 AT1G01020.1 AT1G01020.2 AT1G01030.1 AT1G01040.1 AT1G01040.2 AT1G01046.1 AT1G01050.1 AT1G01060.1 AT1G01060.2 $ cat exon_list.txt | cut -f1 | sort | uniq > gene_name.list
遺伝子は、以下のように40745種類存在することがわかります。
$ cat gene_name.list | wc -l 40745
次に、Genesというフォルダを作り、
$ mkdir Genes
以下のようなコマンドを、
$ cat exon_list.txt | grep AT1G01010.1 > Genes/AT1G01010.1.txt
awkを用いて遺伝子名に対応するように作成し、
“makeGenes.sh”という名前で保存します。
$ head Gene_name.list | awk '{print "cat exon_list.txt | grep "$1" > Gene/"$1".txt"}'
cat exon_list.txt | grep AT1G01010.1 > Genes/AT1G01010.1.txt
cat exon_list.txt | grep AT1G01020.1 > Genes/AT1G01020.1.txt
cat exon_list.txt | grep AT1G01020.2 > Genes/AT1G01020.2.txt
cat exon_list.txt | grep AT1G01030.1 > Genes/AT1G01030.1.txt
cat exon_list.txt | grep AT1G01040.1 > Genes/AT1G01040.1.txt
cat exon_list.txt | grep AT1G01040.2 > Genes/AT1G01040.2.txt
cat exon_list.txt | grep AT1G01046.1 > Genes/AT1G01046.1.txt
cat exon_list.txt | grep AT1G01050.1 > Genes/AT1G01050.1.txt
cat exon_list.txt | grep AT1G01060.1 > Genes/AT1G01060.1.txt
cat exon_list.txt | grep AT1G01060.2 > Genes/AT1G01060.2.txt
$ cat transcript_name.list | awk '{print "cat exon_list.txt | grep "$1" > Genes/"$1".txt"}' > makeTrans.sh
以下のように40745行のシェルスクリプトができるので、
$ wc -l makeGenes.sh 40745 makeGenes.sh
それを実行します。
$ sh makeGenes.sh
すると、以下のように遺伝子ごとに分かれた40745個のファイルができます。
$ ls Genes | head AT1G01010.1.txt AT1G01020.1.txt AT1G01020.2.txt AT1G01030.1.txt AT1G01040.1.txt AT1G01040.2.txt AT1G01046.1.txt AT1G01050.1.txt AT1G01060.1.txt AT1G01060.2.txt $ ls Genes | wc -l 40745 $ cat Genes/AT1G01010.1.txt AT1G01010.1 Chr1 3631 3913 AT1G01010.1 Chr1 3996 4276 AT1G01010.1 Chr1 4486 4605 AT1G01010.1 Chr1 4706 5095 AT1G01010.1 Chr1 5174 5326 AT1G01010.1 Chr1 5439 5899
上記の3列目と4列目の数値群の最小値と最大値が、遺伝子AT1G01010.1の
開始位置、終了位置ということになります。これを統計解析ソフトRで
対話的に行うと以下のようにして求められます。
$ cd Genes/ $ R > in_f <- "AT1G01010.1.txt" > out_f <- "AT1G01010.1.out.txt" > data <- read.table(in_f, sep="\t", quote="") > pos <- c(data[,3], data[,4]) > write.table(cbind(data[1,1:2],min(pos),max(pos)), out_f, sep="\t", append=F, quote=F, row.names=F) > q() Save workspace image? [y/n/c]: y $ cat AT1G01010.1.out.txt V1 V2 min(pos) max(pos) AT1G01010.1 Chr1 3631 5899
これを入力ファイル名と出力ファイル名を指定して、
Rのスクリプトで書くと以下のようになります。
$ vi GeneExt.R args <- commandArgs(trailingOnly = T) in_f <- args[1] out_f <- args[2] data <- read.table(in_f, sep="\t", quote="") pos <- c(data[,3], data[,4]) write.table(cbind(data[1,1:2],min(pos),max(pos)), out_f, sep="\t", append=F, quote=F, row.names=F) q()
GeneExt.Rを以下のコマンドを用いて実行します。
$ R --vanilla --slave --args AT1G01020.1.txt AT1G01020.1.out.txt < GeneExt.R
すると以下のような結果が得られます。
$ cat AT1G01020.1.txt AT1G01020.1 Chr1 8571 8737 AT1G01020.1 Chr1 8417 8464 AT1G01020.1 Chr1 8236 8325 AT1G01020.1 Chr1 7942 7987 AT1G01020.1 Chr1 7762 7835 AT1G01020.1 Chr1 7564 7649 AT1G01020.1 Chr1 7384 7450 AT1G01020.1 Chr1 7157 7232 AT1G01020.1 Chr1 6437 7069 AT1G01020.1 Chr1 5928 6263 $ cat AT1G01020.1.out.txt V1 V2 min(pos) max(pos) AT1G01020.1 Chr1 5928 8737
ここで、awkを使って以下のようなスクリプトを書き、
“GeneExt.sh”という名前で保存します。
$ cat gene_name.list | awk '{print "R --vanilla --slave --args "$1".txt "$1".out.txt < GeneExt.R"}' | head
R --vanilla --slave --args AT1G01010.1.txt AT1G01010.1.out.txt < GeneExt.R
R --vanilla --slave --args AT1G01020.1.txt AT1G01020.1.out.txt < GeneExt.R
R --vanilla --slave --args AT1G01020.2.txt AT1G01020.2.out.txt < GeneExt.R
R --vanilla --slave --args AT1G01030.1.txt AT1G01030.1.out.txt < GeneExt.R
R --vanilla --slave --args AT1G01040.1.txt AT1G01040.1.out.txt < GeneExt.R
R --vanilla --slave --args AT1G01040.2.txt AT1G01040.2.out.txt < GeneExt.R
R --vanilla --slave --args AT1G01046.1.txt AT1G01046.1.out.txt < GeneExt.R
R --vanilla --slave --args AT1G01050.1.txt AT1G01050.1.out.txt < GeneExt.R
R --vanilla --slave --args AT1G01060.1.txt AT1G01060.1.out.txt < GeneExt.R
R --vanilla --slave --args AT1G01060.2.txt AT1G01060.2.out.txt < GeneExt.R
$ cat gene_name.list | awk '{print "R --vanilla --slave --args "$1".txt "$1".out.txt < GeneExt.R"}' > GeneExt.sh
GeneExt.shを、ファイルGenes内に移動させ、実行します。
$ mv GeneExt.sh Genes/ $ cd Genes/ $ sh GeneExt.sh
終了後、コマンドcatを用いて“.out.txt”を含むファイルをすべて
マージし、コマンドgrepを用いて”AT”で始まっている行だけを集めて、
Genes.txtを得ます。
$ cat *.out.txt | grep “^AT” > Genes.txt
以下のように、各遺伝子の染色体上の位置情報のリストが得られました。
$ head Genes.txt AT1G01010.1 Chr1 3631 5899 AT1G01020.1 Chr1 5928 8737 AT1G01020.2 Chr1 6790 8737 AT1G01030.1 Chr1 11649 13714 AT1G01040.1 Chr1 23146 31227 AT1G01040.2 Chr1 23416 31120 AT1G01046.1 Chr1 28500 28706 AT1G01050.1 Chr1 31170 33153 AT1G01060.1 Chr1 33666 37840 AT1G01060.2 Chr1 33666 37780
(4) ゲノムのDNA塩基配列から遺伝子のDNA塩基配列を切り出し、
マルチfasta形式のファイルを作る
最後に、 Genes.txtの遺伝子の位置情報をもとにID_seq.txtから、各遺伝子配列データを切り出し、マルチ fasta形式のファイルを作るスクリプトを作成し、実行します。これから後の操作は、原稿の枚数が多くなってきたので、これは読者への宿題としておきま しょう。
ヒントとしては、PerlやRubyによるスクリプトを用いて作成することも可能ですし（Perlですとsubstr関数を用いてできますし、Rubyですとstr[2..4]のように添字を用いて簡単にできます。）、cutなどのシェルのコマンドをオプション“-c”用いて作成することも可能です。これら一連の操作は、当然のことながら、シェルスクリプトやPerlなどの単一スクリプト言語を用いてすべて記述し本格的なプログラムとすることも可能です。実際に、いろいろな言語で同様の動作をするスクリプトを作成し、速度を比較したりして最適化することもあります。
しかし、ここで紹介した簡便な方法でいろいろな言語のスクリプトをシェルスクリプトで組み合わせることにより、思いついた操作をすぐにスクリプトで書けます。したがって、計算量が少なくすぐに答えを出したい場合には、重宝します。
まとめ
以上のように、いろいろなプログラミング言語や解析ツールをシームレスに組み合わせて、情報を抽出することが可能です。これらの積み重ねは、本格的なビッグデータ解析にもつながっていき、効率的な解析システムの確立にも生かされるものと考えています。
———————————————————————————————————————-
■東京農工大学農学系ゲノム科学人材育成プログラム
OSC2013 Tokyo/Springに出展します。
ぜひお立ち寄りください。
▼セミナー詳細
http://www.ospn.jp/osc2013-spring/modules/eguide/event.php?eid=37
▼展示一覧
http://www.ospn.jp/osc2013-spring/modules/article/article.php?articleid=1