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の必須オプションは以下の通りです。
- INPUT=File
- OUTPUT=File
- RGLB=String
- RGPL=String
- RGPU=String
- RGSM=String
Picardに比べるとだいぶ手間と時間がかかりますが、samtools mergeで追記する方法もあります(参照ページ)。
ご参考までに、手順は以下の通りです。
$ perl -e 'print "\@RG\tID:sample\tSM:sample\tLB:sample\tPL:Illumina\n"' > rg.txt
$ cat rg.txt
@RG ID:sample SM:sample LB:sample PL:Illumina
$ samtools merge -rh rg.txt output.bam sample.bam
【追記】 samtoolsにreheaderという機能が追加されたので、merge機能を使わなくても簡単にread groupを置換できるようになりました。