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

vcf fileをGT(genotype)の情報のみを取り出す

vcf fileからGT(genotype)の情報のみを取り出し、

hom_RR -> 0
het_RA -> 1
hom_AA -> 2

の値に変換するためのスクリプト
あらかじめ、biallelicにフィルタリングしてあります。

vcftools_0.1.13

### GT のみを抽出

vcftools　--gzvcf　test.vcf.gz --extract-FORMAT-info GT --out test

test.GT.FORMATというファイルが作成されます。

このファイルにはGT(genotype)情報のみが記載されています。

scoreに変換して、tsv fileで出力

sed "s/1|1/2/g" test.GT.FORMAT | sed "s/0|1/1/g" | sed "s/1|0/1/g" | sed "s/0|0/0/g" > test_score.tsv

sed "s/1[/]1/B/g" test.GT.FORMAT | sed "s/0[/]1/H/g" | sed "s/1[/]0/H/g" | sed "s/0[/]0/A/g" | sed "s/[.][/][.]/NA/g" > test_score.tsv

間違いなどがありましたら、お知らせいただければ助かります。