kuroの覚え書き

96の個人的覚え書き

vcfにアノテーションを付ける

さて、無事samtools mpileup | bcftools callにてsnv callingが実施できたところで、これを解析するためのアノテーションをつけようと思う。
今回は自前で用意したcontrolのサンプルで検出されたsnp/indelをバックグラウンドのSNVとして対象のvcfにマーキングするのと、ゲノムポジションからarabidopsisのGeneIDに紐付けするためTAIR10のbedファイルを参照することにする。

いろいろ試行錯誤をしながら行きつ戻りつしているが、controlサンプルのSNVをアノテーションに使うためにvcf-annotateを使うことにする。vcf-annotateはVCFToolsの一部分なのでまずはこれをインストール。brewでいける。
そして
controlサンプルのvcfファイルからbedファイルを作成する。
これはawkを使ってvcfファイルから必要な情報を抜いてきて整形し直すことで作成する。

$ cat control.filter.vcf.txt | awk '{OFS="\t"; if (!/^#/){print $1,$2-1,$2,$1":"$2}}' > control.bed

例えばこんな感じ。SNPに正規の名前がないので、仮にchr:posをSNP_IDとしておく
これを圧縮した上で、indexを作っておく。

$ bgzip control.bed
$ tabix -p bed control.bed.gz

ここまで用意できたら、実際にvcfファイルにアノテーションを付ける。

vcf-annotate "$output_dir"/"$samp"_"$exp"_filter.vcf.txt -a control.bed.gz -c CHROM,FROM,TO,INFO/SNP_ID -d key=INFO,ID=SNP_ID,Number=1,Type=String,Description='SNP site' > "$output_dir"/"$samp"_"$exp"_filter_snp.vcf.txt

そして、次にGene_ID をつけるのもvcf-annotateで行う。ここでTAIR10のアノテーションbedファイルを使用すればいいのだけれど、実はCHROMの書式がそのままではchr1,chr2,chr3,chr4,chr5,chrM,chrPとなっていて、今回作成したvcfではそこが1,2,3,4,5,Mt,Ptなので整合性が取れずエラーが出た。
なのでまずはsedで書式を置換しておく。

sed -e "s/chr1/1/g" TAIR10.bed > TAIR10m.bed

こんな感じにすべてのCHROMを書き換えておく。その上で

vcf-annotate "$output_dir"/"$samp"_"$exp"_filter_snp.vcf.txt -a ref/arabidopsis_thaliana/TAIR10m.bed.gz -c CHROM,FROM,TO,INFO/GENE -d key=INFO,ID=GENE,Number=1,Type=String,Description='Gene_ID' > "$output_dir"/"$samp"_"$exp"_filter_snp_anno.vcf.txt

とやればOK

これらの操作によってvcfファイルのINFOカラムに新たなID(SNP_ID, GENE)が追加された。

こんな感じにどんどんとアノテーションを継ぎ足していくことが可能なわけだ。こういう操作はデータをexcelに落としちゃうと手間がかかって全然はかどらないのでvcfファイルの状態でやるに限るね。なので、他に人がやった解析結果をexcelで渡されても本当に困っちゃう。結局最初から解析をやり直すはめになるのよね。<今ここw