以下は個人的なメモ(覚え書き)になります。
GATKでは、BAM fileにRead Group (@RG)が付いてないとエラーが出て解析ができないようです。でもRead Groupって何、どうやって付けるの? 何を書けばいいの? という点で疑問だったので、ちょっと勉強してみました。
まず、以下のページの日本語での解説によると、
http://blog.amelieff.jp/?eid=226752
BAM fileのheaderへのRead groupの入れ方としては、BWAの-R optionを使用するようです。またPicardを使って入れることもできます。実際には、両方やってみました。
ref=/PATH/to/human_g1k_v37_decoy.fasta
reads1=/PATH/to/QC.1.trimmed.fastq
reads2=/PATH/to/QC.2.trimmed.fastq
sample=C8W82ANXX_PG1144_02A16_H1_L003
bwa mem -t 12 -M \
-R "@RG\tID:FLOWCELLID\tSM:${sample}\tPL:illumina\tLB:${sample}_library_1" \
${ref} ${reads1} ${reads2} | samtools view -bS - > ${sample}.bam
ここで”C8W82ANXX_PG1144_02A16_H1_L003"というのは今回の自分のsample file (FASTQ)の名前で、特に何でもいいようです。とりあえずsample fileの名前にしています。
sample=C8W82ANXX_PG1144_02A16_H1_L003
java -jar picard.jar AddOrReplaceReadGroups \
INPUT=${sample}.sort.bam \
OUTPUT=${sample}.sort.addRG.bam \
RGID=FLOWCELLID \
RGLB= ${sample}_library1 \
RGPU=H0164ALXX140820.2 \
RGPL=illumina \
RGSM=${sample}
VALIDATION_STRINGENCY=LENIENT \
CREATE_INDEX=true
Read group (@RG)をこのように付けることで、GATKはエラーを返さずに動くようになりました。
では、Read groupの各項目は何を意味しているのでしょうか?
Read groupの各項目の説明は、SAM/BAM fileについて説明したPDF(samtoolsのHPから)に書かれています。以下、GATKを使うときに入力が求められる5項目についての説明です。
項目 |
説明 |
ID |
SAM/BAM fileの識別名(identifier, ID)。各@RG行は、固有のIDを付けなければならない。 |
LB |
ライブラリ(何だろう? ライブラリの通し番号? 上記のblogの例では"sample"と、別の本の例では"library1"など) |
PU |
Plathome unit(Illuminaの場合はflowcell-barcode-lane, Soildの場合はslide) |
PL |
配列を読むときに使用したプラットフォーム。有効(valid)な値はCAPILLARY, LS454, ILLUMINA, SOLID, HELICOS, IONTORRENT, ONT, PACBIOとのこと |
SM |
サンプルの名前 |
項目PL以外は、特にこの値でなければいけない、ということもなさそうなので、適当に入れれば良い?