2016年2月8日月曜日

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


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