GATKのエラー(2)

GATKでいつも問題になってしまうのが、read group fieldsです。

詳しい説明は以下のページを参照していただければと思います。
https://software.broadinstitute.org/gatk/guide/article?id=6472

BAM fileからのread groupの確認方法は

samtools view -H sample.bam | grep '@RG'　　　　　　　　　　　　

GATKでの解析に必要なread group field は

ID = Read group identifier 
この情報はuniqueである必要があります。
IDが同じの場合、多サンプルでSNP callしても、GATKでは1サンプル分の結果しか出ないという問題が生じます。

PU = Platform Unit 

{FLOWCELL_BARCODE}.{LANE}.{SAMPLE_BARCODE}. 

SM = Sample

PL = Platform/technology used to produce the read 

PL=illumina

と入力しています。

LB = DNA preparation library identifier 

これらの情報はPicardのAddOrReplaceReadGroups を使って記載します。

http://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups

GATKのエラー(1)

GATK実行中に以下のようなエラーが出ました。

SAM/BAM/CRAM file <filename> appears to be using the wrong encoding for quality scores

解決方法は以下のページに書かれています。
https://software.broadinstitute.org/gatk/guide/article?id=6470

古いデータセットでは

Q0 == ASCII 33

ではなく、

Q0 == ASCII 64

となっていることが原因です。
quality scoresが間違って認識されています。

Trimmomaticでも
Paired End Mode:

java -jar <path to trimmomatic.jar> PE [-threads <threads] [-phred33 | -phred64] [-trimlog <logFile>] <input 1> <input 2> <paired output 1> <unpaired output 1> <paired output 2> <unpaired output 2> <step 1> ...

コマンドのうち、[-phred33 | -phred64]を指定するところも間違っていたということになります。

そこでTrimmomaticもやり直しました。

GATKの実行時に
 --fix_misencoded_quality_scores / -fixMisencodedQuals
のオプションを使う必要があります。
https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_engine_CommandLineGATK.php#--fix_misencoded_quality_scores

【170308 追記】
fastq file の時点で、quality scoresを修正しておくことが最も簡単は方法でした。
最初は自分でscoreを33 -> 64に修正していたのですが、DDBJからではなく、ENAからダウンロードすれば既にscoreが修正されていることに気づきました。

quality scoresを確認する最も簡単な方法はfastqcで解析することです。

例えばあるイネのfastq fileですが、
ENAからダウンロードしたファイルでは

Illumina 1.9 ASCII 33になっており、通常の解析で問題ありません。

一方DDBJからダウンロードしたサンプルは

Illumina 1.5 ASCII 64になっており、変換が必要です。

同じ公共データーベースからのダウンロードでも、ファイルが異なることがあります。

確認する必要があるかもしれません。