「オープンソース」を使ってみよう (第27回 フリーソフトウェアを用いたゲノム科学におけるビッグデータ処理)

No Comments
このエントリーをはてなブックマークに追加

—————————————————————————-
東京農工大学ゲノム科学人材育成プログラム特任教授
石井 一夫
—————————————————————————-

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など
 ペタバイト級、エクサバイト級のスケールのビッグデータ処理も視野に入れ
 た次世代ファイルシステムの導入も進んでいます。私たちの教育プログラム
 でゲノムデータ解析に使用している主なフリーソフトを表1に示します。

表:東京農工大学ゲノム科学人材育成プログラムでデータ解析に用いるソフトウェア群
Table:Software for Genomic Data Analysis used in Tokyo University of Agriculture and Technology

 用途  ソフトウェア
 OS  Linux/UNIX
 (CentOS 6.3,
  Scientific Linux 6.3 and FreeBSD 9)
 プログラミング言語  Perl, Python, Ruby,
 Java, C, C++
 データベース  MySQL, PostgreSQL
 ゲノム配列データのアセンブリ  Velvet, ABySS, SOAPdenovo,
 WGS Assembler, MIRA3, Phrap
 ゲノム配列データのマッピング  Bowtie, Bowtie2, BWA, SOAP
 RNA 発現解析用ソフト  Tophat, Cufflinks, Trinity,
 Oases, SOAPdenovo-Trans
 ChIP-Seq解析用ソフト  MACS, Quest, SISSRs, SPP
 統計解析ソフト  R/Bioconductor, Octave
 相同性解析、注釈付けソフト  BLAST, BLAT
 総合DNA解析ソフトウェア  EMBOSS
 生物学データ解析用ライブラリ  BioPerl, BioRuby, BioPython, BioJava
 ビッグデータ解析用ソフトウェア  Hadoop, OpenStack, Eucalyptus

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 <シロイヌナズナの第2番染色体ゲノムの塩基配列>
Chr3 <シロイヌナズナの第3番染色体ゲノムの塩基配列>
Chr4 <シロイヌナズナの第4番染色体ゲノムの塩基配列>
Chr5 <シロイヌナズナの第5番染色体ゲノムの塩基配列>
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

Comments are closed.