さて、無事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で渡されても本当に困っちゃう。結局最初から解析をやり直すはめになるのよね。＜今ここｗ
