Palmsonntagmorgen

NGSデータを使った解析と、その周辺。

マッピング: CRAM形式を試す

Genome Linux NGS

マップデータの形式にはSAM, BAMの他にCRAMという形式があります。

https://www.ga4gh.org/news/cram-compression-for-genomics/

CRAMはBAMと比べて更に高圧縮率だそうです。
今まであまり使う機会が無かったのですが、Twitterで Ewan Birneyさんが強く推していたので、今日はCRAMを試してみようと思います。

SAMtoolsのインストール

今はSAMtools でCRAMも扱えるんですね。便利ですね。
SAMtoolsのインストールは↓を参照してください。

SAMtoolsとリダイレクト - Palmsonntagmorgen

手元のSAMtoolsのバージョンは1.9です。

$ samtools --version
samtools 1.9
Using htslib 1.9
Copyright (C) 2018 Genome Research Ltd.

ゲノムデータの準備

CRAMの変換にはゲノム配列データが必要になります。ゲノム配列作成は以下を参照してください。

常染色体と性染色体のみのゲノム配列ファイル genome.fa を作成する - Palmsonntagmorgen

また、ゲノムに indexが付加されていることが推奨なので、faidxコマンドでindexを作成しておきます。

$ samtools faidx genome_hg19.fa

genome_hg19.fa.fai が作成されます。

bamのダウンロード

bamファイルは何でもいいのですが、今回は以下からダウンロードした GSM733643_hg19_wgEncodeBroadHistoneK562Pol2bStdAlnRep2.bam を使います。

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM733643

（リンク先下部の"Supplementary file" からBAMがダウンロードできます）

BAM -> CRAMの変換

SAMtoolsを使ってBAMからCRAMへ変換してみます。上でindexを作成しましたが、-Tで指定するのはゲノム配列であることに注意してください。 ここではtime コマンドを追加して実行時間を計測しています。
bamのファイル名が長いと見づらいので、ここでは仮にPol2.bamとします。

$ time samtools view -C Pol2.bam -T genome_hg19.fa > Pol2.cram

real    1m46.909s
user    1m37.133s
sys 0m5.343s

変換に時間がかかるという噂がありましたが、シングルコアでも１分ほどで終わりました。

$ ls -lh
合計 906M
-rwxrwxrwx 1 rnakato rnakato 618M  4月 10 17:13 Pol2.bam
-rwxrwxrwx 1 rnakato rnakato 288M  4月 10 17:15 Pol2.cram

確かに、ファイルサイズが半分以下になっています。情報量が同じでサイズが半分以下なのは良いですね。

CRAM -> BAMの変換

ファイルが可逆変換か確認するために、CRAMからBAMに変換してみます。 元のBAMファイルと同一のファイルが生成されることを期待しています。

$ time samtools view -b Pol2.cram > Pol2.bam2
real    1m19.342s
user    1m17.821s
sys 0m1.135s

BAMに戻す時にはゲノムを指定する必要はないようです。

$ ls -lh
合計 1.6G
-rwxrwxrwx 1 rnakato rnakato 618M  4月 10 17:13 Pol2.bam
-rwxrwxrwx 1 rnakato rnakato 683M  4月 10 17:17 Pol2.bam2
-rwxrwxrwx 1 rnakato rnakato 288M  4月 10 17:15 Pol2.cram

何故だか元bamファイルよりファイルサイズが大きくなっています。

３つともsamファイルに変換して調べたところ、Pol2.bam2はPol2.bamと比べてカラムが2つ増えていました。詳細をきちんと確認していませんが、このあたりは元SAM・BAMファイルをどのツールで生成したかによって結果が変わるかもしれません。 なお、Pol2.bam2とPol2.cramは同一でした。

CRAMのソート

$ samtools sort Pol2.cram -O cram -@ 12 > Pol2.sort.cram

問題なくソートできました。 -O cramで出力をCRAMに指定するのを忘れるとbamになってしまうので注意。

bowtie2からのリダイレクト

最後に、bowtie2でマッピングした結果をリダイレクトして直接ソート済CRAMを生成することにチャレンジ。

まず元fastqのダウンロード（何でもよいです）

$ fastq-dump SRR227359

bowtie2でマッピングしてみます。bowtie2のデフォルト出力はSAM形式です。

# SAM
$ bowtie2 -x bowtie2-indexes/UCSC-hg19 -p 24 SRR227359.fastq > SRR227359.sam
# sorted BAM
$ bowtie2 -x bowtie2-indexes/UCSC-hg19 -p 24 SRR227359.fastq | samtools sort > SRR227359.sort.bam
# sorted CRAM
$ bowtie2 -x bowtie2-indexes/UCSC-hg19 -p 24 SRR227359.fastq | samtools sort -O cram - -T genome_hg19.fa > SRR227359.sort.cram
[bam_sort_core] merging from 5 files and 1 in-memory blocks...
[E::cram_get_ref] Failed to populate reference for id 0
samtools sort: failed to write header to "-"

CRAMはsamtools sortに直接投げるとエラーになりました。 修正して以下のように２段階のリダイレクトにします。

# sorted CRAM（修正）
$ bowtie2 -x bowtie2-indexes/UCSC-hg19 -p 24 SRR227359.fastq \
  | samtools view -C - -T genome_hg19.fa \
  | samtools sort -O cram \
  >  SRR227359.sort.cram

今度はうまく行きました。

BAMの時と比較しても実行時間はそれほど変わりません。計算時間的な意味でのオーバーヘッドはほとんど気になりません。

考察

SAMtoolsを使うとCRAMはほとんどストレスなく取り扱えます。ゲノム配列を指定するのがちょっと面倒なくらいです。 計算時間はほとんど変わらず、ファイルサイズは小さくなるので、有用性は高いです。
あとは各種解析ツールが入力にCRAM形式をどの程度受け付けるかによりますね。今のところは対応しているツールはそれほど多くないと思うので、「しばらく使う予定は無いが消すことはできない」というような大量のBAMファイルがある場合にCRAM形式にしておくと容量が削減できてうれしい、というところでしょうか。

あとはpaired-endの時にどうなるのかなどチェックした方がいいかもしれません。