アメリエフのブログ

バイオインフォマティクスの紹介と社員の日々
<< 物事の関係 | main | UniAnno 販売開始キャンペーン! >>
bamにread groupを追記する
GATKは、BAMのフォーマットに厳しく(参照ページ)、たとえばヘッダにサンプル名を含むread groupのリストがあり、かつすべてのリードがそのread groupに属しているBAMしか受け付けません。

Read group(以下RG)は、たとえばBWAではマッピングのときに -R で指定することができます。
$ bwa mem -M -R "@RG¥tID:sample¥tSM:sample¥tPL:Illumina" genome.fa sample_R1.fastq.gz sample_R2.fastq.gz

マッピングの時にRGを付けていなくても、Picard toolkitのAddOrReplaceReadGroupsというツールで追記することができます。

inputのsample.bamにRGを追記してoutput.bamを作成する場合の例を見てみましょう。
$ java -jar AddOrReplaceReadGroups.jar I=sample.bam O=output.bam RGID=sample RGLB=sample RGPL=Illumina RGPU=run_barcode RGSM=sample

AddOrReplaceReadGroups.jarの必須オプションは以下の通りです。

1. INPUT=File
2. OUTPUT=File
3. RGLB=String
4. RGPL=String
5. RGPU=String
6. RGSM=String

Picardに比べるとだいぶ手間と時間がかかりますが、samtools mergeで追記する方法もあります(参照ページ)。
ご参考までに、手順は以下の通りです。
inputのsample.bamにRGを追記してoutput.bamを作成する場合の例を見てみましょう。

まず、以下のようなタブ区切りのテキストファイルを作成します。このファイルの内容がBAMのRGヘッダになります。
$ cat rg.txt
@RG sample SM:sample LB:sample PL:Illumina
rg.txtはこんなコマンドで簡単に作成できます。
$ perl -e 'print "¥@RG¥tID:sample¥tSM:sample¥tLB:sample¥tPL:Illumina¥n"' > rg.txt

次に、samtools mergeには複数のbamファイルが必要なので、ヘッダだけのtmp.bamを用意します。
$ samtools -Hb sample.bam > tmp.bam

最後に、samtools mergeで、sample.bamとtmp.bamをマージさせつつ、rg.txtから新しいbamのRGヘッダを付けます。
$ samtools merge -rh rg.txt output.bam sample.bam tmp.bam
アライメント中のRGタグは、samtools mergeのinputファイル名から付けられます(ここではsample.bamからsampleになります)。このタグがrg.txtで付けたヘッダと一致していないとGATKに怒られるので注意が必要です。

Picardの必須オプションにあまり使いどころのない項目がいくつかありますが、それを省略したいときにこの方法が使えると思います(たとえばここでの例のようにPUを省略してもGATKは動きます)。
| kubo | 次世代シーケンサー解析 | 15:33 | comments(0) | - |









   1234
567891011
12131415161718
19202122232425
262728293031 
<< March 2017 >>

このページの先頭へ